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Preface 


This  report  Is  the  result  of  a  twelve  week  study  on  the  feasibility 
of  using  the  method  of  weighted  residuals  to  determine  approximations  to 
the  discrete  Green's  function  or  an  analog  to  It.  Included  In  the 
report  are  derivations  of  the  methods  of  Galerkin,  collocation  and 
finite  differences,  for  the  one-  aad  two-dimensional  Poisson's  equation. 
The  analytical  solutions  for  various  inhomogenelty  terms  are  also 
included.  The  problem  of  ill-conditioned  matricies  which  arose  in  two 
cases  is  discussed  in  Appendix  A.  All  of  the  primary  goals  of  the  study 
were  met  or  explained. 

During  the  twelve  week  period,  I  learned  a  great  deal  about  the 
theory  of  Green's  functions,  numerical  methods,  matricies  and  the 
problems  that  can  come  about  from  solving  matrix  equations. 

I  would  like  to  acknowledge  Dr.  Kaplan  for  his  support,  and 
seemingly  never-ending  list  of  reference  sources.  His  suggestions 
always  opened  up  new  avenues  for  research,  and  taught  me  more  about  the 
subjects  than  I  really  wanted  to  know. 

I  must  also  thank  my  wife,  Roxann,  for  her  understanding  during 
this  stressful  period,  and  for  being  willing  to  share  the  computer  with 
me  despite  needing  it  for  work  oa  her  own  thesis. 


Dean  E.  Oyler 
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Abstract 


The  purpose  of  this  *tudy'  was  to  determine  the  feasibility  of  using 
the  method  of  weighted  residuals  to  obtain  approximations  to  the 
discrete  Green's  function,  or  analogs  to  it.  The  weighted  residual 
methods  of  Galerlcin  and  collocation,  as  well  as  the  finite  difference 
method  were  programed  on  a  Kaypro  II  micro-computer  in  Microsoft  Basic. 
These  programs  were  used  to  generate  approximations  to  the  one-  and  two- 
dimensional  Poisson's  equation.  The  two-dimensional  case  was  restricted 
to  the  geometry  of  a  unit  square.  Various  inhomogeneity  terms  were  used 
ro  obtain  approximate  solutions  to  the  discrete  Green's  functions  or 
their  analogs.  The  results  were  compared  with  the  analytical  values  at 
various  interior  nodal  points  on  the  mesh.  The  average  percent  error 
for  the  approximations  were  reported  for  each  case  as  the  number  of 
interior  nodal  points  was  Increased.  The  areas  of  consideration  were 
the  rate  of  convergence  to  the  analytical  solution,  the  amount  of  time 
it  took  to  run  each  program,  and  the  accuracy  of  the  approximate 
solutions.  The  results  of  this  study  indicate  that  the  Green's 
functions  or  analogs  obtained  are  valid  approximations  to  the  discrete 
Green's  functions.  The  method  of  weighted  residuals  proved  to  be  very 
sensitive  to  the  choice  of  basis  functions,  resulting  in  ill-conditioned 
matricles  in  some  Instances. 


INVESTIGATION  OF  THE  NUMERICAL  METHOD  OF  MOMENTS  FOR 
DIGITAL  COMPUTER  DETERMINATION  OF  DIFFERENTIAL  EQUATIONS 


I.  Introduction 


Background 


The  final  step  in  the  mathematical  treatment  of  many 
ph}3ics  and  engineering  is  often  finding  the  solution  to 
value  problem.  The  standard  differential  equations 
encountered  in  mathematical  physics  include  - 


problems  in 
a  boundary- 
most  often 


LaPlace's  equation: 

V2#  *  0 


(1) 


Poisson's  equation: 

V24>  -  -  p 


(2) 


the  wave  equation: 

V2*  -  i/c2[afy/*t2]  -  o  (3) 


and  the  Helmholtz  equation: 

(V2  +  k2)*  -  f  U) 

Frequently,  the  solutions  to  these  equations  can  be  represented  in 
terms  of  a  Green's  function.  There  are  several  advantages  in  the  use  of 
Green's  functions  as  solutions  to  these  boundary-value  problems.  One 
advantage  is  that  it  enables  a  differential  equation  with  suitable 


boundary  conditions  to  be  solved  by  an  ordinary  integral.  Another 
advantage  is  that  once  the  Green's  function  for  a  particular 
differencial  operator  and  geometry  has  been  found,  it  can  be  utilized 
for  all  other  problems  involving  the  same  differential  operator  and 
geometry,  but  with  different  expressions  for  the  inhomogeneity  or  source 
term.  If  these  Green's  functions  for  different  differential  operators 
and  geometries  could  be  tabulated,  they  could  be  used  to  solve  boundary-* 
value  problems  quite  easily,  in  a  manner  analogous  to  the  use  of  a  table 
of  Integrals. 

Although  Green's  functions  have  been  obtained  analytically  for 
certain  standard  geometries  (planes,  rectangles,  spheres,  cylinders), 
for  the  usual  equations  of  mathematical  physics  (£qs(l-4)),  a  difficulty 
arises  in  finding  the  Green's  functions  for  mixed  or  Irregular 
geometries.  In  these  situations,  one  must  employ  the  use  of  numerical 
methods  techniques. 

Previous  thesis  research  and  publications  have  successfully  solved 
the  problem  of  numerically  approximating  Green's  functions  by  means  of 
finite  difference  algorithms  (1,3,6).  In  this  method,  one  uses 
approximations  of  derivatives  (usually  a  truncated  Taylor's  series)  to 
convert  the  boundary-value  problem  into  a  large  series  of  simultaneous 
algebraic  equations,  which  can  than  be  solved  with  relative  ease  using 
matrix  methods  on  a  digital  computer. 

This  thesis  will  investigate  the  use  of  a  different  numerical 
technique,  the  method  of  weighted  residuals,  to  solve  the  necessary 
differential  equations.  In  this  technique,  the  unknown  solution  is 
expressed  as  a  series  of  functions  which  can  be  manipulated  to  once 
again  reduce  the  problem  to  solving  a  series  of  simultaneous  algebraic 


equations.  The  inverse  of  the  coefficient  matrix  of  these  equations  is 
analogous  to  the  discrete  Green's  function  for  the  differential 
operator. 


Objective 

This  research  effort  will  compare  the  discrete  Green's  functions, 
obtained  by  both  the  method  of  weighted  residuals  and  the  method  of 
finite  differences,  for  both  the  one-  and  two-dimensional  case,  with 
ra°oect  to  accuracy  and  feasibility  for  digital  computation. 

The  current  thesis  problem  is  a  follow-up  to  a  previous  M. S.  thesis 
(1)  that  reported  conflicting  results  for  the  one-  and  two-dimensional 
cases,  concerning  which  of  the  three  methods  was  best.  This  study  will 
attempt  to  verify  or  refute  the  conflict  between  the  two  cases  by 
recreating  parts  of  the  previous  thesis  using  a  different  computer  coda, 
and  different  matrix  solving  routines. 

Scope 


This  study  will  only  consider  tne  problem  of  the  one-  and  two- 
dimensional  Poisson's  equation,  with  Dirichlet  boundary  conditions.  Tha 
solutions  for  both  the  one-  and  two-dimensional  Green's  functions  will 
be  compared  using  the  finite  difference  method  and  the  method  of 
weighted  residuals. 

Approach 


The  initial  approach  to  this  3tudy  will  be  to  develop  computer 
programs  which  use  the  method  of  finite  differences  and  two  of  the 


methods  of  weighted  residuals  (Galerlcin's  method,  aad  collocation)  to 
obtain  approximations  for  both  the  discrete  Green's  function,  or  its 
analog,  and  the  solution  for  the  one-dimensional  Poisson's  equation. 
Homogeneous  Dirichlet  boundary  conditions  will  be  assumed  for  all  cases. 

Once  the  initial  programs  have  been  developed,  they  will  then  be 
modified  to  handle  the  two-dimensional  cases. 

The  usefulness  of  the  Green's  functions  obtained  in  the  previous 
steps  will  then  be  analyzed  by  varying  the  inhomogeneity  term  for  the 
Poisson's  equation.  Areas  of  consideration  will  include  the  number  of 
calculations  required,  computer  run  time  for  each  method,  the 
convergence  rate  to  the  correct  solution,  the  overall  accuracy  of  the 
approximations  as  compared  to  the  analytical  solution,  and  how  the 
results  compare  to  the  earlier  study  (1). 

Finally,  the  feasibility  and  possible  directions  for  continued 
research  into  these  approximation  methods  will  be  explored. 


POISSON'S  EQUATION  IN  ONE-DIMENSION 


The  initial  problem  examined  in  this  study  is  the  one-dimensional 
Poisson's  equation.  The  general  form  of  the  problem  can  be  expressed  as 


L  u(x)  -  g(x) 


2  2 

where  L  is  the  linear  differential  operator,  d  / dx  ,  g(x),  the 

Inhomogeneity  term,  is  the  source  or  excitation  (a  known  function),  aad 
u(x)  is  the  field  or  response  (the  unknown  function  to  be  determined) 
(7:1-2). 

Associated  wifch  the  problem  in  this  study  are  the  homogeneous 
Dlrichlet  boundary  conditions 


u(0)  «  0 
u( 1 )  -  0 


(6-a) 

(6-b) 


Analytical  Solution 


The  general  solution  to  Eq(5),  with  3n  inhomogeaei ty  term  of  the 


g(x)  ■  Ax*  +  Bx  +  C 


can  be  found  by  direct  Integration  to  be 


u(x)  «  Ax^/12  +  3x3/6  +  Cx2/2  +  Dx  +•  £ 


By  applying  the  boundary  conditions  (Sq  (6)),  Eq  (3)  becomes 


u(x)  -  Ax^/12  +  Bx3/6  *  Cx2/2  -  (A+2B*SC)x/12 


(9) 


which  is  Che  analytical  solution  to  the  one-dimensional  Poisson's 
equation. 

Numerical  Approximations 

All  of  the  numerical  approximations  in  this  study  make  use  of  a 

technique  in  which  a  mesh  is  superimposed  over  the  region  of  interest  of 

the  problem.  For  the  one-dimensional  case,  this  merely  Involves 
subdividing  the  region  by  M  equally  spaced  interior  nodes.  The 
numerical  method  is  then  applied  at  these  nodes  resulting  in  a  set  of  N 
simultaneous  algebraic  equations,  which  can  be  solved  to  give  the 

approximate  solution  to  the  problem. 

The  accuracy  of  the  approximation  depends  upon  the  number  of 

interior  nodal  points  used.  A  fine  mesh  with  many  nodal  points  will 
generally  result  in  a  more  accurate  solution  to  the  problem.  A  sample 
mesh  for  t he  one-dimensional  case  is  shown  below  in  Figure  1 


Figure  1.  1-0  Mesh  With  5  Interior  Modes 
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The  mesh  in  Pigure  1  has  five  interior  nodes,  with  a  step  size  h  between 
each  of  the  nodes  equal  to  1/(N+1),  to  ensure  that  they  are  equally 
spaced. 

Plni te  Difference  Method.  The  method  of  finite  differences  makes 
use  of  a  truncated  Taylor's  series  to  approximate  the  solution  to  the 
problem.  According  to  Taylor's  theorem  (14:6),  when  a  function  u,  and 
its  derivatives  are  singla-valued,  finite  and  continuous  functions  of  x, 
then 

u(x+h)  -  u(x)  +  hu'(x)  +  h2a"(x)/2  +  hJu'"(x)/6  +  ...  (10) 

and 

u(x-h)  -  u(x)  -  hu'(x)  +  h2u"(x)/2  -  bV"(x)/6  +  ...  (11) 

Adding  these  two  expressions  yields 

2 

u(x+h)  +  u(x-h)  »  2u(x)  +  h  u"(x)  +  higher  order  terms  (12) 

Neglecting  the  higher  order  terms  and  solving  for  u"(x) 

u"(x)  »  d2u/dx2  -  (u(x+h)  -  2u(x)  +  u(x-h))/h2  (13) 

This  equation  allows  an  approximation  to  Eq  (5)  to  be  made  so  that 

(u(x+h)  -  2u(x)  +  u(x-h))/h2  -  g(x)  (14) 

By  applying  this  approximation  to  each  nodal  point  in  the  mesh,  the 
equation  becomes  a  series  of  N  simultaneous  algebraic  equations  which 
can  then  be  solved  by  matrix  techniques.  These  equations  can  be  written 


in  matrix  notation  as 


C  u  -  g  (15) 

In  all  cases  when  vector  notation  is  used,  a  boldfaced  capital  letter 
indicates  a  rectangular  matrix,  and  a  boldfaced  lower  case  letter 
indicates  a  column  vector. 

The  coefficient  matrix,  C,  formed  from  these  equations  is  known  as 
a  tri-diagonal  matrix.  A  tri-diagonal  matrix  has  non-zero  values  only 
along  the  main  diagonal,  and  the  adjacent  diagonals  both  above  and  below 
the  main  diagonal.  As  an  example,  the  tri-diagonal  matrix  equation  for 
five  interior  nodes  is  shown  in  Figure  2. 
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Figure  2.  Tri-diagonal  F.D.  Matrix  Equation  for  5  Modes 


Method  of  Weighted  Residuals.  In  the  various  methods  of  weighted 
residuals,  the  unknown  function  u(x)  from  £q(5)  is  expanded  in  a  series 
of  functions,  b^(x),  b2(x),  b^(x),  ...  in  the  domain  of  L,  as 

u(x)  b  (*)  (16) 

T  n  Q 

where  the  a  are  constants,  and  the  b  (x)  are  expansion  or  basis 
n  n 

functions  (7:5-6).  The  basis  functions  are  chosen  so  as  to  match  the 
boundary  conditions  of  the  problem.  For  the  boundary  conditions  of 


8 


Eq(6),  the  basis  functions  were  chosen  to  be  a  power  series  (7:7)  of  the 


b  (x)  •  x  -  x 


n  *  1 ,2 ,3 ,  .  .  . 


By  substituting  Eq(16)  into  £q(5)  and  due  to  the  linearity  of  L,  the 
equation  can  be  rewritten  as 

2»aL  bQ(x)  -  g(x)  (18) 

ft 

Sow,  a  set  of  weighting  or  testing  functions  w^ ,  »  w3>*  *  *  defined 

in  the  range  of  L  (7:6),  and  the  inner  product  of  Eq(13)  and  the 

weighting  functions  w  is  taken  so  that 

ni 


a  <w  ,  L  b  (x)>  -  <w  ,  g(x)> 
nan  a 


or  in  the  more  condensed  matrix  notation 


m  *  1,2,3,  «  •  • 


C  i  *  g 


where  G  is  the  square  coefficient  matrix 


’<w1,  L  bj(x)>  <w^  L  b2(x)>  .  . 
<»2,  1  Oj(x)>  <w2*  L  b2(x)>  »  •  • 


a  is  the  column  vector 


and  g  13  che  column  vector 


<w1.  g1<x» 
<»2i  S2U>> 


The  column  vector  of  constants,  a,  can  be  determined  by 


a  -  C-1  g 


Once  a  has  been  determined,  the  approximate  solution  can  then  be  written 


u(x)  *  b*  a 


for  a  method  valid  over  the  entire  region,  or  as 


a  «  B  a 


for  a  method  valid  only  at  discrete  points,  where  b*  is  the  transpose 


vec  tor 


[bl»  b2’  b3»  *•*» 


a(x1) 


u(x2} 


.u(VJ 


^2^X1  ^ 9  •••* 


b,  (x^),  b0(xM) ,  «•••  bM(x%J)| 


The  particular  choice  made  for  the  weighting  functions  w^ 
determines  which  of  the  methods  of  weighted  residuals  Is  being  used. 

Galerkln^s  Method.  The  choice  of  letting  the  weighting 

3H”1 

function  w_  ■  b  (x)  Is  known  as  Galerkin's  method.  For  w  «  x  -  x  , 
m  n  m 

the  values  for  the  coefficient  matrix,  C  can  be  found  by  taking  the 

Inner  product  of  w_  and  L  b  (x).  The  result  Is 

in  u 

l 

C  «  <w  ,  L  b(x)>  -  f  (x  -  X**1)  d2/dx2(x  -  x1**1)  dx  (30) 
Similarly,  g  Is  found  to  be 

®m  *  g(*)>  *  J "(*  “  x°+1)  d2/dx2(Ax2  +  Bx  +  C)  dx  (31) 

• 

There  are  no  limitations  on  the  value  of  x,  therefore  the  results 
obtained  by  Galerkin's  method  are  valid  over  the  entire  region,  and  not 
just  at  the  nodal  points  of  the  mesh. 

Collocation  Me thod .  This  procedure,  also  known  as  point- 

matching  (7:10),  makes  use  of  the  Dirac  delta  function  as  the  weighting 

function.  The  Inner  product  of  w  and  L  b  (x)  Is  then 

m  n 

C  -  <w  ,  L  b  (x)>  »  f  6(x  -  x  )  d2/dx2(x  -  xn+1)  dx  (32) 
<nn  m  n  j  m 

e 

where  x,  are  the  points  equispaced  In  the  Interval  0  <  x  <  1, 

ul 

(x^  *  m/(M+l),  m-1,2,3,  .  .  S)  corresponding  to  the  nodal  points.  And 

gm  Is  given  by 

*m  *  <wm»  S(*)>  *  / 8(x  -  x^,)  d2/dx2(Ax2  Bx  ¥  C)  dx  (33) 

o 

The  collocation  method  limits  the  value  of  x  to  the  values  at  the 


nodal  points,  therefore,  the  method  of  collocation  approximations  are 
only  valid  at  these  nodal  points. 

Greenes  Pune tions  and  Analogs 

This  saction  will  show  how  the  various  approximation  techniques 
utilized  in  this  study  may  be  related  to  the  Green's  function  or  its 
analog. 

Finite  Difference  Method.  The  relationship  for  the  Green's 
function  for  the  differential  operator  L  of  Eq C 5 )  has  been  defined 
(15:6-7)  to  be 

d2G(x|x  )/dx2  -  8(x-x  )  (34) 

0  o 

where  G(x|xq)  is  the  Green's  function  for  Eq(5),  x  is  the  field  point, 
and  xq  is  the  source  point.  The  associated  Dirichlat  boundary  conditions 
for  Eq(34)  are 


G(0lx  )  -  0  (35-a) 

0 

G(l|x0)  -  0  (35-b) 

When  applied  to  a  mesh  with  step  size  h^,  Eq(34)  takes  the  form  of 

the  discrete  Green's  function  (5:314-315) 

d2G  (x|x  )/dx2  -  8(x-x  )/h  (36) 

NO  ON 

where  N  is  once  again  the  number  of  interior  nodal  points  on  the  mesh. 

The  derivative  term  can  then  be  treated  in  the  same  manner  as  was 
done  earlier  for  the  finite  difference  method,  and  be  replaced  by  a 
central  difference  equation 


i 


l 


i 


i 


I 


I 


i 


« 


« 


12 


(37) 


d2SN(x|xQ)/dx2 

-  (GJ1(x+hfJ|xo)  -  2Gn(x|xq)  +  GM(x-h^|xo))/h2 

Substituting  this  equation  into  Eq(36)  and  multiplying  through  by  h4", 
the  expression  becomes 

GS<X+hM!xo)  "  2GH(xlXo}  +  G^(x"hHlXo)  "  ^  8(x-xo)  (38) 

Applying  Eq(38)  along  with  the  associated  boundary  conditions  to 

2 

each  of  the  nodal  points  yields  a  series  of  N  simultaneous  algebraic 
equations  which  can  be  expressed  in  matrix  notation  a3 

0  %  -  “*  <39> 

where  C  is  once  again  the  coefficient  matrix,  G  is  the  discrete  Green's 
function  matrix,  and  1^  is  the  identity  matrix  of  order  N.  fhe 
coefficient  matrix  of  the  finite  difference  method  of  Eq(15)  is 
equivalent  to  the  inverse  of  the  numerical  Green's  function  matrix 
multiplied  by  h.  The  finite  difference  equations  and  -the  Green's 
function  equations  yield  Identical  approximations  to  the  solutions. 

He thod  of  Weighted  Residuals.  For  the  one-dimensional  Poisson's 
equation  (Eq(5>),  and  its  associated  boundary  conditions  (Sq(6)>,  the 
Green's  function  for  the  problem  can  be  determined  analytically  (15:1- 
12).  The  solution  to  Eq(5)  with  its  various  inhomogeaei ty  terms  can  be 
found  by  calculating  the  integral 

i 

u(x)  «  / G(x|xq)  g(x)  dxQ  (40) 

o 

where  G(x|xQ)  is  the  Green's  function  for  Sq(5)  and  its  associated 


boundary  conditions,  x  is  the  field  point,  and  x  is  the  source  point. 


ID 


Eq(40)  can  be  written  in  matrix  notation,  for  the  discrete  Green' 3 
function  on  a  mesh  of  step  size  h.^  as 

I  •  h  G  t  (41) 

where  a  is  the  column  vector  of  solutions  at  discrete  points  on  the  mesh 
for  the  given  inhomogeneity.  The  tilde  has  been  placed  over  the  left- 
hand  side  of  the  equation  to  stress  the  fact  that  the  discrete  Green' s 
function  solution  may  not  necessarily  be  equal  to  the  solution  obtained 
by  the  method  of  weighted  residuals  for  a. 

If  Eq( 24)  is  substituted  into  Eqs(25)  and  (26),  they  can  be  written 
as 

u  -  b*  C"1  g  (42) 

a 

and 

Uy  «  B  C  1  g  (43) 

These  two  equations  for  the  weighted  residual  approximations  do  not 
contain  the  factor  for  the  step  size,  n^,  as  does  Eq(41)  This  is  due  to 
the  fact  that  the  weighted  residual  approximations  involve  a  summation 
over  tJ  terms,  while  the  discrete  Green's  function  solution  was  developed 
by  approximating  an  Integral  equation  (Eq(4G)),  where  the  step  size,  h^, 
corresponds  to  the  dx  term.  The  analog  to  the  discrete  Green's  function 
in  Eqs(42)  and  (43)  can  be  defined  to  be 

G*  -  b*  C~l  (44) 

for  a  method  valid  over  the  entire  region  of  interest,  such  as  Galer- 
kin's  method  (the  bar  over  the  G  is  used  to  indicate  a  column  vector  in 
this  case  to  avoid  confusion  with  the  inhomoganei ty  terra),  and 
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*  -1 

G,  -  B  C  (4 

for  a  method  valid  only  at  discrete  points,  such  as  the  collocation 
method.  The  superscript  asterisk  Indicates  that  they  are  analogs  to  the 
discrete  Green's  function.  They  are  considered  analogs  since  the 
elements  that  make  up  the  Inhomogeneity  vector,  g,  In  Eq(23)  may  not 
necessarily  be  equal  to  the  inhomogeneity  term  g(x)  because  of  the 
weighting  factor. 

Computer  Analysis 

All  of  the  numerical  approximation  routines  were  developed  on  a 
Kaypro  II  microcomputer  in  Microsoft  Basic,  using  double-precision 
arithmetic.  The  program  listings  are  included  in  Appendix  C. 

Inhomogeneity  Terms  and  Exact  Solutions.  The  same  four  functions 
chosen  by  Clapp  (1)  for  the  excitation  terms  were  also  used  in  this 
study  so  that  comparisons  could  be  made  with  the  results  obtained  in 
that  earlier  study.  These  four  functions  were 


1)  L  u(x)  -  10  (46-a) 

2)  L  u(x)  -  x1 2 3 4  (46-b) 

3)  L  u(x)  -  x2  +  1  (46-c) 

4)  L  u(x)  ■  x2  +  x  +  1  (46-d) 

with  the  associated  boundary  conditions 

u(0)  -  0  (47-a> 


The  number  of  interior  nodes  were  carefully  chosen  so  that  the 
solution  comparisons  could  all  be  made  at  the  same  nodal  points,  x»l/3 


and  x»2/3. 

The  analytical  solutions  to  the  problem  set  were  found  by  direct 
integration  to  be 


1) 

u(x) 

-  5x2  - 

5x 

(43-a) 

2) 

u(x) 

-  x*/12 

-x/12 

(43-b) 

3) 

u(x) 

-  x4/12 

-  x2/2  -  7x/12 

(43-c) 

4) 

u(x) 

-  x4/12 

+  x^/3  +  x2/2  -  3x/4 

(43-d) 

The  exact  solutions  at  the  comparison  points  are  listed  in  Table  1 


TABLE  1 

Exact  Solutions  to  Eq(43) 

Problem  #  1  x  »  1/3  x  «  2/3 


-1.111111  -1.111111 

-.026749  -.039095 

-.137360  -.150260 

-.137243  -.211934 


Average  Error.  The  average  percent  error  was  the  criterion  by 
which  the  correctness  of  the  approximations  were  measured.  For  the  one¬ 
dimensional  case,  the  average  percent  error  was  defined  to  be 


<EkI>  -  |uVI(l/3)  -  u(  1/3)|  +  I uVT ( 2/ 3 )  -  u( 2/3)1  )  100 

s  -8~  ~Wi) -  ~8~h/3) -  }•  — 


where  u^  is  the  approximation  at  a  specific  point  for  a  given  number  of 
nodes,  and  u  is  the  exact  solution. 


Comparison  of  Approxlma tions  to  Exact  SolutloQ.  The  plots  of  the 
average  percent  error  vs.  number  of  interior  nodal  points  are  shown  in 
Figures  3-6  for  each  of  the  four  equations  in  the  problem  sat.  fhe 
actual  values  of  the  approximations  are  included  in  Appendix  B. 

Except  in  two  instances,  both  of  the  methods  of  weighted  residuals 
(Galerkin  and  collocation)  yielded  approximations  to  the  exact  solution 
that  were  orders  of  magnitude  better  than  the  method  of  finite 
differences.  In  addition  they  were  able  to  achieve  these  good 
approximations  using  a  relatively  low  number  of  interior  nodal  points, 
whereas  the  finite  difference  approximations  converged  to  the  exact 
solution  more  slowly.  The  only  deviations  from  thi3  trend  were  in  the 
case  of  problem  1  (Eq(46-a))  when  all  methods  did  equally  well  on  the 
average,  and  in  problem  4  (Eq(46-d)),  when  approximations  in  the 
Galerkin  routine  began  to  diverge  rapidly  from  the  correct  solution  for 
eight  or  more  interior  nodes. 

Compu ter  Run  Times.  Each  program  was  timed  by  hand  using  an  elec¬ 
tronic  stop-watch  to  obtain  values  for  the  computer  run  time.  The  times 
varied  by  only  one  or  two  seconds  for  the  various  inhomogeneity  terms, 
therefore  the  times  reported  are  averages  for  each  method.  The  plot  of 
program  run  times  V3  number  of  interior  nodes  is  shown  in  Figure  7. 

In  all  cases,  the  method  of  finite  differences  took,  the  least 
amount  of  time  to  run.  This  method  not  only  has  fewer  intermediate 
calculations  than  the  other  methods,  but  in  addition,  the  coefficient 
matrix  was  tri-diagonal;  this  allowed  the  use  of  an  extremely  efficient 
coutine  expressly  written  to  30lve  tri-diagonal  matricies  (2:122). 

3oth  methods  of  weighted  residuals  took  considerably  longer  to 
arrive  at  a  solution  (in  some  cases,  10  to  30  times  longer),  for  a  given 
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Average  Error  (;*>) 


Figure  3.  Average  Error  vs  dumber  of  Interior  nodes 
for  g(x)  •  10 
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Figure  4.  Average  Erro^  va  Humber  of  Interior  Sode3 
for  g(x)  -  x 
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Figure  5.  Average  Error  vs  Number  of  Interior  Nodes 
for  g(x)  -  x^  +  1 
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Average  Error  (tf) 


Legend 

O  —  Pinite  Difference 
Q  — Galetkln 

A —  Collocation 


11  14  17 


Number  of  Interior  Nodes 


Figure  6.  \verage  Error  vs  dumber  of  Interior  Nodes 
for  g(x)  «  x  t  1 


number  of  nodal  points.  The  collocation  program  was  the  slowest  of  the 
three  methods,  largely  due  to  the  matrix  solving  technique  U3ad  In  this 
program.  This  program  was  based  on  an  Inversion  routine  using  Causs- 
Jordan  elimination  with  column  shifting  (3:294-295),  whereas  the 
Galerkin  program  required  the  use  (see  Appendix  A)  of  a  direct  Saussian 
elimination  routine  with  partial  pivoting  (9:192-193).  Inversion 
methods  require  more  computations  to  arrive  at  the  solution  than  do 
direct  methods,  hence  the  large  difference  In  the  times  between  the  two 
methods  of  weighted  residuals. 

Overall  Solution  Accuracy,  and  Comparison  with  Earlier  Results. 
The  results  of  this  study  do  not  support  those  reported  in  the  earlier 
study  (1).  Clapp  reported  that  In  all  cases,  the  finite  difference 
method  was  superior  to  the  method  of  weighted  residuals  for  the  one- 
JQ<  dimensional  case,  yielding  average  errors  of  about  one  percent  when  17 
or  more  nodes  were  used.  This  study  shows  that  both  methods  of  weighted 
residuals  achieved  average  errors  on  the  order  of  10  ^  percent  after 
only  five  interior  nodes  were  used.  This  accuracy  was  matched  by  the 
finite  difference  method  only  for  a  constant  Inhomogeneity  term. 

done  of  the  oscillations  reported  by  Clapp  for  the  collocation 
method  were  noted  in  this  study. 

Conclusions 

From  the  analysis  performed,  it  is  clear  that  the  method  of 
weighted  residuals  is  the  best  choice  if  few  Interior  nodes  are  desired. 
While  these  methods  are  more  difficult  to  program  and  take  longer  to  run 
than  the  method  of  finite  differences,  the  overall  accuracy  is  much 
superior. 
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III.  POISSOcTS  EQUATION  Id  TWO  DIMENSIONS 

The  other  problem  examined  In  this  study  is  the  two-dimensional 
Poisson's  equation.  The  general  form  of  the  problem  can  be  expressed  in 
the  same  manner  as  Eq(5) 


L  u(x,y)  -  g(x,y)  (50) 

2  2 

where  L  is  now  the  two-dimensional  linear  differential  operator,  d  /dx 
2  2 

"*■  d  / dy  ,  g(x,y)  is  the  two-dimensional  inhomogenei ty  term,  and  u(x,y) 
is  the  unknown  function  to  be  determined.  The  Dirichlet  boundary 
conditions  associated  with  the  two-dimensional  problem  are 


u(0,y) 

-  0 

(51-a) 

u(l,y) 

-  0 

( 51-b) 

u(x,0) 

-  0 

(51-c) 

u(x,l) 

-  0 

( 51-d) 

Equations  (50)  and  (51)  define  Poisson's  equation  for  the  region  of  a 
unit  square. 


Analytical  Solution 


The  analytical  solution  to  Eq  (50),  with  its  associated  boundary 
conditions  Eq  (51),  can  be  found  by  a  Fourier  series  expansion  (6:41-42) 
to  be 


u(x,y) 


►  /W2  2  2  /  /  dtj  sin(mirx)sin(niry)3in(ni7r£  )sin(nirT|  )g(£,n) 

m-1  n-1  J  J  T - : - 


,  2  .  2v 

(a  +  n  ) 


(52) 


For  a  given  excitation,  g(x,y),  the  series  solution  to  Sq  (52)  can 
be  found  by  integration.  This  solution  can  then  be  programed  on  a 


computer  to  provide  a  numerical  solution,  at  a  specific  point,  for  the 


desired  number  of  summation  terms.  It  should  be  noted  that  because  the 
solution  involves  a  double  summation  process,  obtaining  the  result  may 
take  a  considerable  length  of  time. 

Numerical  Approximations 

For  the  two-dimensional  case,  the  mesh  is  superimposed  over  a  unit 
square,  with  s  equally  spaced  nodes  in  the  x  direction,  of  step  size  b, 
and  t  equally  spaced  nodes  in  the  y  direction,  of  step  size  k.  The 
total  number  of  interior  nodes  in  the  mesh,  N,  is  equal  to  s  X  t.  V 
sample  mesh  for  the  two  dimensional  unit  square,  with  the  interior  nodal 
points  numbered,  is  shown  in  Figure  8. 


Sines  the  number  of  interior  nodal  points  used  determines  the  number  of 
simultaneous  algebraic  equations  that  must  be  solved,  the  x  and  y  step 
sizes  were  chosen  not  to  be  equal.  This  was  done  in  order  to  keep  the 
matricies  down  to  a  manageable  size. 

Fini te  Difference  Method.  As  was  shown  in  section  II  (the  one- 
dimensional  Poisson's  equation)  the  finite  difference  matricies,  and  tha 
Green's  function  matricies  were  identical,  therefore,  the  development  of 
the  finite  difference  method  will  only  be  done  utilizing  the  Green's 
function  approach. 

The  Green's  function  for  the  two  dimensional  linear  operator,  L,  of 
Eq  (50),  is  defined  to  be  (15:1-18) 

i2G(x|xo;y|yQ)/*x2  +  lZG(x|xo;y| yQ)/Jdy2  -  5(x-xo)fi(y-yo)  (53) 
with  the  associated  boundary  conditions  (Eq(51))  being  represented  by 


G(0|xo;y|yo) 

-  0 

(54-a) 

5<1lXo;yiyo) 

-  0 

(54-b) 

G(xl  V°lyo> 

-  0 

( 54-c ) 

G(xlxo;1fyo) 

-  0 

( 54-d) 

The  discrete  Green's  function  imposed  on  a  mesh  with  step  3lze  h  in 
the  x  direction,  and  k  in  the  y  direction,  can  be  written  as  (4:315) 

a23N(x!xo;y|yo)/»x2  +  *23N(x|xo;y|yo)/^y2  -  5(x-xQ)5(y-yo)/hk  (55) 

The  derivative  terms  can  once  again  be  replaced  by  central 
difference  expressions 
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>2V*lv^0)/a*2 

*  (G5l(x+h|xo;y|yo)  -  2GN<x|  xQ ;y|  yQ)  +  G^Cx-hj xq; y| yQ) )/h2  (56-a) 

and 

a2GN(xlxo;y!yo>/,y2 

-  (3N(xjx0;y+kjy0)  -  2G^(x| xQ;y| yQ)  +  GN(x|xQ;y-lt| jq) )/k2  (56-b) 

aad  thea  substituted  back  into  Eq(55),  so  that  the  final  form  of  the 
equation  becomes 

k2G^(x+h|xo;y|yo)  -  2k2GM<x|xo;y|yQ)  +  k2Gs(x-h|xo;y|yo) 

+  h2GJ1(xlXo;y+'<lyo>  "  2h2GN(xlx0;ylyo)  *  h'Gi^xlxo;y"xlyo)  (57) 

-  nk  5(x-xo)8(y-yo) 

Applying  this  equation,  and  its  associated  boundary  conditions 

( Eq( 54) )  at  each  of  the  tf  interior  nodes  of  the  mesh  will  result  in  a 
2 

series  of  N  simultaneous  algebraic  equations.  These  equations  can  be 
represented  in  matrix  notation  as 

C%  9  hkIN  (58) 

where  C  is  the  coefficient  matrix,  is  the  discrete  Green's  function 
matrix,  and  1^  is  the  identity  matrix  of  order  N.  Eq(58)  can  be  solved 
using  matrix  techniques  to  find  G^,  the  approximation  to  the  discrete 
Green's  function  matrix. 

Me thod  of  Weighted  Residuals.  For  the  two-dimensional  method  of 
weighted  residuals,  the  form  of  the  solution  is  almost  Identical  to  the 
one  dimensional  case  (Eq(16)) 

u(x,y)  “2®  b  (x,y)  (59) 

n  n  a 
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where  the  aQ  are  again  constants,  and  the  on(x,y)  are  the  basis 
functions,  which  were  chosen  to  satisfy  the  boundary  conditions 

it 

(Eq(51)).  Ozisik  (13:340-344)  developed  basis  functions  for  a 

rectangular  region  as  the  product  of  a  function  f(x,y),  and  various 
powers  of  x  and  y.  For  the  geometry  of  the  problem  in  this  study,  the 
function  f(x,y)  was  chosen  to  be 


f(x,y)  -  xy(l-x)(l-y) 


and  the  basis  functions  were  chosen  to  be  the  same  as  those  used  by 


Clapp  (1) 


b  »  (xy)2^a  f (x, y)  (for  n-1,4,7,...) 

n 

-  x2(n+l ) /3  (for  a«2,5 ,8 , . . . ) 

bQ  «  y2a/3  f (x,y)  (for  n*2,8,9,...) 


(61-a) 

(61-b) 

(61-c) 


The  basis  functions  were  split  into  these  three  subgroups  ia  order 
to  facilitate  the  necessary  integrations  that  follow  in  the  sections  oa 
the  Galerkin  and  collocation  methods. 

By  substituting  Sq(59)  into  Eq(5Q),  the  equation  becomes 

2*  Lb  (x,y)  *  g(x,y)  (°2) 

n 

As  in  the  one-dimensional  case,  a  set  of  weighting  functions,  w^, 
is  then  defined  in  the  range  of  L,  and  the  inner  product  of  Eq(62)  and 
these  weighting  functions  is  taken  30  that 


Lbn(x,y)>  -  <w!n,  3(x,y)>  m-1,2,3,... 


which  C3n  again  be  represented  in  the  more  condensed  matrix  notation  as 


where  C  is  Che  square  coefficient  matrix 


<wlf  Lb1(x,y)>  <wx,  Lb2(x,y)>  .  .  . 
<w2,  Lbj (x,y)>  <w2,  Lb2(x,y)>  .  .  . 


(65) 


a  is  the  column  vector 


(66) 


and  g  is  the  column  vector 


<wx.  g1(*»y)> 

<w2»  S2<‘X,y)> 


(67) 


The  values  of  a  can  then  be  computed  from 


a 


g 


and  the  approximate  solutions  can  then  be  written  as 

u(x,y)  -  bf  a 


(68) 


(69) 


for  a  method  valid  over  the  entire  region,  or  as 


for  a  method  valid  oaly  at  discreta  points,  where 


b  ^  bi *  b2’  b3’  *  *  *  * 

u(x1,y1) 

u(x2,y2) 


Lu(vVJ 


b1(*l,y1),  b2(x1,y1),  o^(xl,y1) 

»  -  ’>2U2’72) . WV  (73) 

bI b2 VXS,y?P 

Galerkin"s  Method.  In  the  two-dimensional  method  of  Galerkin, 
the  weighting  functions  are  again  chosen  as  being  equal  to  the  basis 
functions  Eq(61).  The  values  for  the  coefficient  matrix,  C  can  be  found 
by  taking  the  inner  product  of  wfl  and  L  ba(x,y),  resulting  in 


Gmn  *  <VV  Lbn^x»y^>  “  Jj  +  ^ bn(x,y)  dx  dy  (74) 


and  g  can  ba  found  in  a  similar  manner  to  be 
m 


8„  *  <w|B,  g(x,y)> 


I  l 

ffo  (d^/dx^  +■  d^/dy^)  (Ax^  +  By^  +■  Cx  +  Dy  +  E)  dx  dy  (75) 


whera  In  each  case,  wm  is  defined  to  be  equal  to  t>n(x,y)  ( Eq (51))  for 
the  various  values  of  n. 

Thera  are  no  restrictions  on  the  value  of  x  or  y  in  2qs(74)  and 
(75),  so  the  Salerkin  approximations  are  valid  over  the  entire  region. 

Colloca  tion  Me  thod.  In  the  two-dimensional  method  of 
collocation,  the  weighting  functions  are  chosen  to  be  equal  to  the  two- 
dimensional  Dirac  delta  function 


w  -  8(x— x  )8(y— y_)  (76) 

in  m  m 

where  the  coordinate  x  is  defined  to  ba  the  x  coordinate  of  the  nth 

m 

interior  node  and  y  is  defined  to  be  the  y  coordinate  of  the  mth 

m 

interior  node.  The  values  for  the  coefficient  matrix,  C  can  be  found  by 
taking  the  Inner  product  of  w^  and  L  bQ(x,y) 


c™  -  <v 

i  i 

-  Hr y„) 

o  o 

and  g  can  be  found  to  be 


(d2/dx2  +  d2/dy2)  b^(x, y)  dx  dy 


(77) 


gm  *  <wm’  g(*’y)> 


S' 

o  « 


8(x-x  )8(7~7_)  (d2/dx2  +  d2/dy2)  (Ax2  ♦  By2  *■  Cx  (73) 

fa  ul 

+  Dy  +■  £)  dx  dy 


The  values  of  x  and  y  ara  restricted  to  the  coordinates  of  the 
interior  nodal  points  for  the  collocation  method,  therefore,  the 
approximations  are  only  valid  at  these  points,  3nd  not  over  the  entire 
region. 
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Green's  Pune tlons  and  Analogs 

Since  the  finite  difference  method  was  developed  in  terms  of  the 
Green's  function,  only  the  method  of  weighted  residuals  will  be 
addressed  in  this  section. 

For  the  two-dimensional  Poisson's  equation  (£q(50)),  and  its 
associated  boundary  conditions  (Eq(51)),  the  Green's  function  for  the 
problem  can  be  determined  analytically  (6:42-43).  The  solution  to 
Eq(50)  with  its  various  inhomogeneity  terms  can  be  found  by  calculating 
the  Integral 

i  I 

u<x.y)  mJJG(x\*0i y  1  y0>  g(*.y>  dxo  dy0  (79) 

o  O 

where  G(x|xo;y|yQ)  is  the  Green's  function  for  £q(50)  and  its  associated 
boundary  conditions,  x  and  y  ara  the  field  point  coordinates,  and  xq  and 
yQ  are  the  source  point  coordinates. 

Eq(79)  may  be  written  in  matrix  notation,  for  the  discrete  Green's 
function  on  a  mesh  of  step  sizes  h  In  the  x  direction,  and  tc*  In  the  y 

3  C 

direction,  as 

%  *  **  GN  *N  (80) 

where  <1  is  the  column  vector  of  solutions  at  discrete  points  on  the  mesh 
for  the  given  inhomogeneity,  and  the  tilde  again  signifies  that  the 
discrete  Green's  function  solution  may  not  necessarily  be  equal  to  the 
solution  obtained  by  the  method  of  weighted  residuals,  u. 

The  two-dimensional  forms  of  Eqs(42)  and  (43)  are  tne  approximate 
solutions 

-  b’  <f 1  g  (31) 
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and 


B  C  g 


( 


Agaia,  the  weighted  residual  approxima tions  do  not  contain  the  factors 
for  the  step  sizes  h  and  k.  This  is  due  to  the  fact  that  the  weighted 
residual  approximations  involve  a  summation  over  N  terms,  while  the 
discrete  Green's  function  solution  was  developed  by  approximating  an 
integral  (Eq(79))t  where  the  step  sizes  h  and  Ic  correspond  to  the  dx  and 
dy  terms.  The  analog  to  the  discrete  Green's  function  in  Eqs(31)  and 
(82 )  <-an  be  defined  to  be 


b  G_1 


(33) 


for  a  method  valid  over  the  entire  region,  such  as  Galerkin's  method 
(where  the  bar  notation  again  indicates  a  column  vector,  and  was  used  to 
prevent  confusion  with  the  inhomogeneity  term),  and 


G*  «  B  <f 1  (84) 

for  a  method  valid  only  at  discrete  points,  such  as  the  collocation 
method.  The  asterisk  indicates  that  they  are  analogs  to  the  discrete 
Green's  function,  and  that  the  elements  that  make  up  the  Inhomogeneity 
vector,  g,  may  not  necessarily  be  equal  to  the  lnhomogeneity  term  g(x,y) 
because  of  the  weighting  factor. 


Computer  Analysis 

For  the  two-dimensional  case,  new  programs  were  developed  from  the 
one-dimensional  programs  to  handle  the  approxima tions  for  the  three 
techniques.  The  program  listings  are  Included  in  Appendix  C. 
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Inhomogeneity  Taras  and  Exact  Solutions.  The  same  four  excitation 
terms  chosen  by  Clapp  (1)  were  used  la  the  two-dimensional  case  so  that 
comparisons  could  be  made  with  the  results  obtained  la  his  study.  The 
four  problems  were 


1)  L 

u(x)  -  10 

(35-a) 

2)  L 

/  x  2 

u(x)  ■  X 

(35-b) 

3)  L 

2  2 
u(x)  -  x  +  y 

(35-c) 

4)  L 

2  2 

u(x)  ■  x  +  y  x 

(35— d) 

wl  th  the 

associated  Oirichlet 

boundary  conditions 

u(0,y)  -  0 

(36-a) 

u(l,y)  -  0 

(36-b) 

u(x,0)  -  0 

(86-c) 

u(x,l)  »  0 

(36-d) 

The 

number  of  interior 

nodes  were  carefully  chosen  so 

that  the 

solution 

comparisons  could 

all  be  made  at  the  same  four 

x,y  nodal 

points;  d/3,1/3),  (2/3, 1/3),  (1/3, 2/3),  and  (2/3, 2/3). 

The  analytical  solutions  to  the  problem  set  were  found  by 
Integrating  Eq(52)  with  a  general  form  of  the  equations  used  for  the 
inhomogenei ty  term, 

g(x,y)  ■  Ax2  +  By2  +  Cx  +  Dy  +  E  (87) 

where  A,  B,  C,  0  and  S  are  all  constants. 

The  solution  for  the  general  Inhomogenei ty  term  with  the  boundary 
conditions  of  Eq(86)  is 
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u(x,y)  -  -4/ir22  2[sin(m,rx>sin(n,ry)/(m2  +  n2)] 
m*l  n»l 

{A[(2/mV  -  l/mir)(-l)in  -  2/mV]  (l/nir)(l  -  (-l)a) 

+  B[(2/n3ir3  -  i/nir)(-l)a  -  2/n3ir3]  (l/mir)(l  -  (-l)"1)  (38) 
-  C((-l)°/aiiwr2l(l  -  (-1)Q)  -  D((-l)a/nniJ^l (1  -  (-l)m) 

+  Etl/aair2l((-l)m  -  l)((-l)a  -D) 

This  solution  was  then  programed  oa  the  computer,  and  the  numerical 
values  for  the  lnhomogenei ties  at  each  point  were  generated  by  summing 
over  both  m  and  n  from  one  to  seventy.  Each  numerical  value  then  was 
obtained  using  4900  summation  terms.  The  numerical  values  for  tha  four 
problems  are  listed  In  Table  2. 


TABLE  2 

Exact  Solutions  to  Eq(38) 


Problem  if 

(1/3, 1/3) 

(2/3, 1/3) 

(1/3, 2/3) 

(2/3, 2/3) 

1 

-.6034615 

-.6034615 

-.6034515 

-.5034515 

2 

-.0126051 

-.0230496 

-.0126051 

-.0230496 

3 

-.0252102 

-.0356547 

-.0356547 

-.0460993 

4 

-.0501611 

-.0710501 

-.0606056 

-.0814946 

Average  Error.  For  the  two-dimensional  case,  tne  average  percent 
error  was  defined  to  be 

<E  >  -  j  [u  .(point  1)  -  u(point  1 )(  (  .  100  (39) 

S  |  1  u(polnt  1)  )  4 

where  u..  Is  the  approximation  at  a  specific  point,  and  u  Is  the  exact 
N 

solution. 
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Comparison  of  Approximations  to  the  Exac t  Solutions.  The  plots  of 
the  average  percent  error  vs  number  of  interior  nodal  points  for  each  of 
the  four  aquations  in  the  problem  set  are  shown  in  Figures  9-12. 

A3  in  the  earlier  study  (1),  tne  collocation  method  failed  to 
converge  to  the  correct  solution  despite  using  a  direct  Gaussian 
elimination  routine  instead  of  performing  a  matrix  inversion.  A  more 
detailed  discussion  of  the  problem  with  the  collocation  method  can  be 
found  in  Appendix  A. 

In  each  of  the  four  cases  for  the  other  two  techniques,  the 
Galerkin  method  yielded  better  results  than  the  finite  difference 
method,  which  is  not  only  consistent  with  the  results  reported  for  the 
one-dimensional  case,  but  also  with  the  results  reported  by  Clapp. 

Computer  Run  Times.  Each  program  was  timed  by  hand  using  an 
electronic  stop-watch  to  obtain  values  for  the  computer  run  time.  The 
times  were  again  averaged,  since  they  differed  by  only  one  or  two 
seconds.  The  plot  of  program  run  time  vs  number  of  Interior  nodal 
points  is  shown  in  Figure  13. 

In  all  cases,  the  method  of  finite  differences  was  again  the 
quickest  of  the  three  programs  to  run.  While  the  two-dimensional  finite 
difference  coefficient  matrix  was  not  a  tri-diagonal  matrix  (as  in  the 
one-dimensional  case),  it  was  still  a  relatively  sparse  one  (ie.  few 
non-zero  terms).  Because  there  were  fewsr  computations  required  to 
create  the  finite  difference  coefficient  matrix  than  for  the  other  two 
methods,  it  ran  faster  (despite  being  based  on  a  matrix  inversion 
routine).  Both  methods  of  weighted  residuals  were  developed  using  the 
direct  Gaussian  elimination  routine. 
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Average  Error  (#) 


Figure  9.  Average  Error  vs  clumber  of  Interior  Nodes 
for  g(x,y)  *  10 
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Figure  10.  Average  Zrror  vs  Number  of  Interior  Nodes 
for  g(x,y)  -  x* 
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Figure  11.  Average  Error  js  Number  of  Interior  tfodes 
for'g(x,y)  -  x  y 


.0 


Average  Error  (%) 


30J 


254 


20  J 


15  4 


10  H 


5  1 


— r- 

10 


Legend 

O-  finite  Difference 
□  -  Galarkin 


— T" 

16 


— T— 

22 


-i — 

28 


Number  of  Interior  Nodes 


— i — 

3^ 


Plgure  12.  Average  Error  vs  riunber  of  Interior  tfodes 
for  g(x,y)  -  x  +  y  +  x 
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Figure  13.  Computer  Aun  Times  vs  dumber  of  Interior  ^odes 
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Overall  Solution  Accuracy,  and  Comparison  with  Earlier  Elesul ts. 
The  results  of  this  section  of  the  study  support  those  reported  by  Clapp 
(1),  the  Galerkin  method  yielded  better  results  than  did  the  finite 
difference  method,  with  initial  average  errors  from  1.5  to  3.0  percent, 
compared  with  13.5  to  27.0  percent.  The  Galerkin  method  converged  to 
within  one  percent  after  34  interior  nodes  were  used,  while  the  finite 
difference  method  remained  above  four  perceat. 

A  slight  oscillation  in  the  Galerkin  results  was  noticed  for  22 
interior  nodal  points. 

Conclusions 

Prom  the  analysis  performed,  the  method  of  weighted  residuals  was 
once  again  superior  to  the  method  of  finite  differences,  for  fewer 
interior  nodal  points. 

0)  *  Meither  method  was  as  accurate  as  its  one-dimensional  counterpart. 

This  is  most  likely  due  to  round  off  error  caused  by  the  larger  matrix 
sizes,  and  to  the  rather  lenythy  recursion  relations  used  to  create  the 
coefficient  matricies. 

While  the  Galerkin  method  wa3  clearly  better  than  the  finite 
difference  method,  it  Ls  also  more  difficult  and  tedious  to  change  the 
inhomogeneity  term,  or  the  form  of  the  differential  equation,  due  to  the 
Integrals  that  must  be  solved.  Despite  this  difficulty,  the  Galerkin 
approach  i3  clearly  the  method  of  choice  for  solving  the  two-dimensional 
Poisson's  equation. 


43 


IV.  Conclusions  and  ilecoamenda  tions 


Conclusions 

Several  points  should  be  made  concerning  the  use  of  the  method  of 
weighted  residuals  for  determining  approximations  to  the  discrete 
Green's  function. 

First,  both  the  Galerkin  and  collocation  methods  yield  analogs  of 
the  Green's  function,  which  are  as  useful  a3  the  discrete  Green's 
function  itself,  and  they  can  be  used  (at  least  theoretically)  to  find 
the  solution  to  the  one-  and  two-dimensional  Poisson's  equation  with 
various  inhomogeneity  terms.  A  major  drawback  to  the  use  of  the  method 
of  weighted  residuals  is  that  the  inhomogeneity  matrix  mu3t  be 
recalculated  for  each  different  inhomogeneity  term.  This  involves  the 
calculation  of  a  lengthy  double  integration  in  the  two-dimensional 
problem,  especially  in  the  case  of  Galerkin' s  method. 

The  next  point  is  the  criticality  of  the  choice  of  basi3  functions. 
Both  the  one-dimensional  Galerkin  routine,  and  the  two-dimensional 
collocation  routine  seemed  to  show  the  effects  of  the  choice  of  basis 
functions  resulting  in  ill-conditioned  matricies. 

Finally,  the  method  of  weighted  residuals  takes  more  time  to  run  on 
the  computer  than  does  the  finite  difference  method.  For  low  numbers  of 
interior  nodes,  this  may  not  be  much  of  a  problem,  but  in  the  two- 
dimensional  programs,  calculations  involving  34  nodal  points  took  25 
minutes  to  arrive  at  a  solution. 

Recommanda  tions 

The  first  step  in  any  follow-on  study  should  be  the  actual 
calculation  of  the  Green's  function  for  the  various  metnods,  so  that 


they  can  be  compared  for  accuracy. 


In  addition,  one  other  major  area  requires  further  study  in 
utilizing  the  method  of  weighted  residuals  -  the  choice  of  basis 
functions.  An  ia-depth  study  of  the  orthogonal  collocation  method,  and 
the  choice  of  basis  functions  in  general  would  be  most  beneficial. 

The  choice  of  using  a  micro-computer  for  developing  the  programs 
was  probably  not  very  wise.  Although  the  results  were  as  accurate  as 
those  done  on  mainframe  computers,  the  programs  took  too  long  to  run. 
Future  work  should  be  done  on  a  larger,  faster  machine. 
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Appendix  A 


Ill-conditioned  Matricles 


The  divergence  from  the  correct  solution  reported  for  the  one- 

2 

dimensional  Galerkin  routine  with  g(x)  ■  x  +  x  +  1,  and  for  the  two- 
dimensional  collocation  routine,  were  both  the  results  of  ill- 
conditioned  matricles. 

Ill-conditioned  matricles  are  nearly  singular  systems  (la.  those 
that  have  no  unique  solution)  which  are  extremely  sensitive  to  small 
changes  in  the  coefficient  matrix,  C,  and  the  right  hand  side,  g.  For 
example,  tne  solutions  to  the  following  two  similar  matrix  equations 
differ  greatly: 


-1  1.00001 


1  -1  *  u j  'l 

-1  1.00001  u,  0 

’ 1  -i  Fuji  r  f 

m 

-1  0.99999  u  0 

»  L  - 


The  solution  to  Eq(90)  is  [10001  100000],  and  the  solution  to  Eq(91)  is 

[-99999  -100000]  (11:343).  This  shows  that  a  relatively  small  change 

in  the  value  of  b7  >00002  can  cause  a  large  change  in  the  solution. 

It  is  easiest  to  Interpret  what  is  happening  geometrically.  Each 
solution  may  be  thought  of  as  representing  the  point  of  intersection  of 
two  nearly  parallel  lines.  Any  slight  shift  in  either  of  the  two  lines 
will  greatly  change  the  point  of  Intersection. 

Ill-conditioning  of  a  system  may  be  attributed  to  any  of  the 


following  sources  (11:343-347): 


1.  For  well-coadi tioaed  physical  problems,  ill-conditioned 
equations  may  be  caused  by  a  correct,  but  very  fine  mesh  idealization. 
Mo  numerical  problems  are  encountered  wnen  the  problem  is  solved  with 
coarse  idealizations.  As  the  mesh  is  repeatedly  subdivided,  the 
condition  number  increases.  Eventually  the  buildup  of  error  due  to 
round  off  in  the  calculations  swamps  any  accuracy  Improvement  due  to  the 
finer  discretization. 

2.  The  form  of  the  right-hand- side  vector  can  have  a  significant 
effect  in  many  applications.  For  example,  in  the  bending  and  stretching 
of  a  flat  plate,  the  stiffness  matrix  may  uncouple  into  an  ill- 
conditioned  submatrlx  and  a  well-conditioned  submatrix.  In  complex 
structural  systems,  weak  coupling  can  occur  so  that  a  force  vector 
acting  on  the  ill-conditioned  part  will  excite  the  ill-conditioning; 
whereas  if  it  acts  on  the  well-conditioned  part,  highly  accurate 
solutions  are  produced. 

3.  The  condition  of  the  system  is  influenced  by  tne  choice  of 
basis  functions. 

The  two-dimensional  collocation  routine  seemed  to  exhibit  the 
traits  listed  in  source  1.  The  program  was  run  for  different  numbers  of 
interior  nodal  points  than  the  ones  used  for  the  other  programs.  Figure 
14  shows  a  plot  of  the  average  percent  error  vs  number  of  Interior  nodal 
points  for  the  collocation  method  when  fewer  nodal  points  were  used. 
The  plot  shows  that  the  collocation  method  was  in  fact  converging  to  the 
correct  solution  until  the  mesh  size  passed  a  critical  value. 

The  one-dimensional  Salerkin  routine  seemed  to  exhibit  the  traits 

of  source  2.  It  worked  perfectly  well  for  all  other  inhomogeneity 

2 

terms,  but  failed  for  g(x)  ■  x  +  x  *■  1  beyond  11  interior  nodal  points. 


Average  Error  (£) 


One  way  of  determining  if  a  matrix  is  ill-conditioned  is  to  observe 

whether  the  condition  number  increases  significantly  as  the  size  of  the 

matrix  is  increased.  The  condition  number  for  a  matrix  can  be 

calculated  by  multiplying  the  value  of  C.  by  C  }  ,  where  C. 

1  f  (&ex  1 }  ttta  x  i  f  iQdx 

is  defined  to  be  the  summation  of  all  of  the  terms  within  the  column  of 

the  matrix  leading  to  the  greatest  value,  and  C~*  is  the  same 

i,max 

calculation  performed  on  the  inverse  of  the  matrix. 

Table  3  shows  how  the  condition  number  for  the  two-dimensional 
collocation  coefficient  matrix  increased  as  the  size  of  the  matrix  was 
increased. 


Table  3 

Condition  lumbers  for  the  2-d  Collocation  Matrix 


Size  of  Matrix 

Condition  Number 

4X4 

2.7  X  101 

5X6 

5.9  X  102 

3X3 

1.6  X  103 

10  X  10 

2.6  X  105 

20  X  20 

12 

3.9  X  10 

Little  can  be  done  to  improve  the  solutions  obtained  from  an  Ill- 
conditioned  matrix.  One  improvement  though  is  to  use  a  direct  method  of 
solving  the  matrix,  rather  than  an  Inversion  routine.  This  was  doae  for 
the  one-dimensional  Galer'xin  program  whan  it  was  discovered  that  the 
approximations  were  all  diverging  from  the  correct  solution  after  only 
eight  Interior  nodes.  Once  the  matrix  inversion  routine  was  replaced 
with  a  direct  Gaussian  elimination  routine,  the  results  Improved 
significantly,  except  in  'the  aforementioned  case. 


Changing  che  matrix  solving  routiae  was  aot  enough  however  to  aalce 


a  difference  In  the  two-dimensional  collocation  program.  Others 
(11:345-347,  12:60-65)  have  suggested  that  the  ill-conditioning  may  be 

the  result  of  the  Improper  choice  of  basis  functions;  that  orthogonal 
polynomials  would  be  better  suited  for  use  in  thi3  situation.  k 
detailed  discrlption  of  the  orthogonal  collocation  method  can  be  found 
in  Plnlayson  (4:97-107). 


Numerical  Approxlma tioa3  for  the  Problem  Sets 


This  Appeadix  contains  the  values  for  the  finite  difference, 
Galerkin,  and  collocation  approximations  to  Eq(46)  and  the  finite 
difference  and  Galerkin  approximations  to  Eq(85)  (the  results  of  the 
two-dimensional  collocation  routine  are  discussed  in  Appendix  A). 

The  values  listed  are  the  outputs  from  the  computer  routines  listed 
in  Appendix  C,  rounded  to  seven  decimal  places.  These  values  are  listed 
by  inhomogeneity  term  within  each  method,  and  are  tabulated  according  to 
the  number  of  interior  nodes  used.  The  average  percent  error  is  also 
listed  for  each  set  of  values. 


Finite  Dlffereace  Approximations  for  g(x)  »  13 


Humber  of 
laterlor 
tfodes 

x  -  1/3 

mm 

Average 

Error 

(Z) 

2 

-1.1111111 

-1.1111111 

0.0 

5 

-1.1111109 

-1.1111109 

l.l  E-5 

8 

-1.1111111 

-1.1111111 

0.0 

11 

-1.1111113 

-1.1111113 

2.1  E-5 

14 

-1.1111110 

-1.1111110 

6.3  E-6 

17 

-1.1111108 

-1.1111103 

2.3  E-5 

Table 

5 

Finite  Difference  Apprgxlmations  for  g(x) 

2 

-  X 

Humber  of 
Interior 

Node  3 

x  -  1/3 

mm 

Average 

Error 

(?) 

2 

-.0246913 

-.0370370 

6.5 

5 

-.0262346 

-.0385802 

1.6 

3 

-.0265203 

-.3386602 

0.7 

11 

-.0266204 

-.0389661 

3.4 

14 

-.0266666 

-.0390123 

0.3 

17 

-.0266913 

-.0390375 

0.2 

Number  of 
Interior 
'lodes 


3 


2/3 


Average 

Error 

CO 


2 

-.1358025 

-.1481431 

1.4 

5 

-.1373457 

-.1496913 

0.4 

8 

-.1376315 

-.1499771 

0.2 

11 

-.1377315 

-.1500772 

0.1 

14 

-.1377777 

-.1501234 

6.0  E-2 

17 

-.1373029 

-.1501486 

4.0  E-2 

Table 

7 

Finite  Difference  Approximations  for  g(x)  »  x' 

'  +  x  4-  1 

dumber  of 
Interior 

Modes 

x  -  1/3 

x  -  2/3 

Average 

Error 

(%) 

2 

-.1351352 

-.2093765 

1.0 

5 

-.1867234 

-.2114197 

0.3 

8 

-.1370142 

-.2117055 

0.1 

11 

-.1371142 

-.2118056 

6.0  E-2 

14 

-.1371605 

-.2118518 

4.0  E-2 

17 

-.1371356 

-.2113769 

■■MB 

Table  3 


Galerkin  Approximations  for  g(x)  ■  10 


Number  of 
Interior 

Nodes 

; 

x  -  1/3 

x  -  2/3 

Average 

Error 

U) 

2 

-i. mini 

-1.1111112 

5 

-l.niini 

-1.1111112 

3 

-l.iiiiin 

-1.1111112 

6.1  E-6 

11 

-l.niini 

-1.1111112 

6.1  E-6 

14 

-l.niini 

-1.1111112 

6.0  E-6 

17 

-l.niini 

-1.1111112 

6.1  2-6 

Table  9 

-  Galerkin  Approximations  for  g(x) 


x 


2 


Number  of 
Interior 

Nodes 

m 

mm 

Average 

Error 

U) 

2 

-.0271605 

-.0395062 

1.3 

5 

-.0267490 

-.0390947 

3.5  E-6 

3 

-.0267490 

-.0390947 

3.5  E-6 

11 

-.0267490 

-.0390947 

3.5  E-6 

14 

-.0267490 

-.0390947 

3.5  E-6 

17 

-.0267490 

-.0390947 

3.5  E-6 

Number  of 
Interior 
Nodes 


3 


2/3 


Average 

Error 

(%) 


2 

-1.1111111 

-1.1111112 

6.1  E-6 

5 

-1.1111111 

-1.1111112 

6.1  E-6 

3 

-1 .1111111 

-1.1111112 

6.1  E-6 

11 

-1.1111111 

-1.1111112 

6.1  E-5 

u 

-1.1111111 

-1.1111112 

6.2  2-6 

17 

-1.1111111 

-1.1111113 

3.3  2-6 

Table 

13 

Collocation  Approximations  for  g(x)  ■  x 

2 

Number  of 
Interior 

Nodes 

x  -  1/3 

mm 

Average 

Error 

(%) 

2 

-.0246914 

-.3703704 

6 . 5 

5 

-.0267490 

-.0390947 

3.5  E-6 

8 

-.0267490 

-.0390947 

3.5  E-6 

11 

-.0267490 

-.0390947 

3.5  E-6 

14 

-.0267490 

-.0390947 

3.5  E-5 

17 

-.0267490 

-.0390947 

6.1  E-6 

dumber  of 
Interior 
Modes 

mm 

■nza 

BBS 

K 

mm 

— W 

Average 

Error 

(*> 

4 

-.5128205 

-.4059829 

-.5341330 

-.5123205 

13.5 

10 

-.5653312 

-.5469123 

-.5634183 

-.5640845 

7.0 

16 

-.5734571 

-.5674212 

-.5744274 

-.5729557 

5.2 

22 

-.5759343 

-.5732613 

-.5763574 

-.5756971 

4.7 

23 

-.5769715 

-.57:,:  534 

-.5771928 

-.5763423 

4.4 

34 

-.5774941 

-.5766637 

-.5776240 

-.5774159 

4.3 

Table  17 


2 

Finite  Difference  Approximations  for  g(x,y)  -  x 


dumber  of 
Interior 
Modes 

R9 

m 

mu 

IBM 

m 

Average 

Error 

(7.) 

4 

-.0092593 

-.0146605 

-.0100309 

-.0135135 

25.3 

10 

-.0114660 

-.0202691 

-.0115141 

-.0210923 

9.4 

16 

-.0118848 

-.0213449 

-.0119352 

-.0216323 

6.1 

22 

-.0120202 

-.0216807 

-.0120430 

-.0213120 

5.1 

28 

-.0120734 

-.0218200 

-.0120906 

-.0213904 

4.6 

34 

-.0121082 

-.0218333 

-.0121155 

-.0219309 

4.4 

d 


€>« 


FINITE  DIFFERENCE  ROUTINE  -  (ONE-DIMENSION) 

10  REM  **^********^*^***********************^ k*lrk*itkk*** »★********»»***■*» 
20  REM  THIS  PROGRAM  UTILIZES  A  ROUTINE  THAT  WAS  TAKEN  FROM  ELEMENTARY 
30  REM  NUMERICAL  ANALYSIS,  BY  CONTE  AND  D£  BOOR,  PG  122;  AND 
40  REM  TRANSLATED  INTO  BASIC.  ALL  VALUES  ARE  IN  DOUBLE- PRECIS ION.  ^ 

60  PRINT" ENTER  N,  THE  NUMBER  OF  STEPS"; 

70  DEFDBL  Q,R,H 
30  INPUT  N 

90  LPRINT"PINITE  DIFFERENCE  ROUTINE  USING  DIRECT,  TRIDIAGONAL  APPROACH" 

100  LPRINT 

110  LPRINT"N«  ";N 

120  Q-N 

130  H«1#/(QH#) 

140  DIM  A#(N),B#(N),C#(N),D#(N) 

150  FOR  I-  1  TO  N 
160  R»I 

170  a:HI)— 1^(Q+1#)"2 
180  C#(I)»A#(I) 

190  D#(I)-2#*(Q+1#)"2 
200  B#(I)»1#+40*U*H)~2 
210  NEXT  I 
220  GOSUB  290 

230  PRINT"TH£  SOLUTION  IS:" 

240  FOR  I-  1  TO  N 
250  PRINT  I,B#(I) 

260  LPRINT  I,B#(I) 

270  NEXT  I 
280  END 

290  IF  N>1  THEN  320 
300  3#(1)-B#(1)/D#(1) 

310  RETURN 

320  FOR  I»  2  TO  H 

330  R— A#(I)/D#(I-1) 

340  D#(I)«  D#( I)+R*C#( 1-1 ) 

350  8*(I)-B#(I)+R*B#(I-1) 

360  NEXT  I 

370  B#(N)»B#(N)/D#(N) 

380  K-N 

390  FOR  J-  2  TO  N 
400  K-K-l 

410  B#(K)-(B#(K)-C#<K)*B#(K-H))/D#(K) 

420  NEXT  J 
430  RETURN 
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FINITE  DIFFERENCE  ROUTINE  -  (TWO-DIMENSION) 


1  REM** ******************************************************** ****** 

2  REM  THIS  PROGRAM  USES  AN  INVERSION  ROUTINE  TAKEN  FROM  NUMERICAL 

3  REM  METHODS  8 f  R.W.  HORNBECK,  PG  294-295.  rH£  ROUTINE  EMPLOYS 

4  REM  GAUSS- JORDAN  ELIMINATION  WITH  COLUMN  SHIFTING  TO  MAXIMIZE 

5  REM  PIVOT  ELEMENTS 

10  DEPDBL  S,T,H,K,A-G,X,Y 
20  PRINT  CHR$(26)+CHR$(27)+CHa$(13) 

30  PRINT-TWO-DIMENSIONAL  FINITE  DIFFERENCE  GREEN'S  FUNCTION  ROUTINE" 
40  PRINT: PRINT SPRINT 

50  PRINT-INPUT  S,  THE  #  OF  X  DIVISIONS;  AND  T,  THE  #  OF  Y  DIVISIONS" 

60  INPUT  S,T 

70  PRINT  S,T 

80  N-S*T 

90  LPRINT"N»";N 

100  H-l#/(S+l#):K»l#/(T+liO 

110  PRINT"H-";H;"K-";K 

120  DIM  G#(N,N),J#(N+25),F#(N),U#(N) 

140  REM  CREATES  INITIAL  COEFFICIENT  MATRIX  FOR  THE  GREEN'S  FUNCTION 
150 

160  FOR  1-1  TO  N 
170  IF  I  MOD  S  -1  THEN  190 
130  G#( 1,1-1 )-K*K 
190  G#(I,I)— 2#*(X*K+H*H) 

200  IF  I  MOD  S-0  THEN  220 
210  G#(I,I+1)-K*K 
220  IF  (l-SX-0  THEN  240 
230  G#(I,I-S)-H*H 
240  IF  ( I+S ) >-N  THEN  260 
250  G#(I,I+S)-H*H 
260  NEXT  I 
270  FOR  1-1  TO  N 
230  FOR  J-l  TO  M 
290  PRINT  G#( I,J) 

300  NEXT  J 
310  PRINT 
320  NEXT  I 

340  REM  MATRIX  INVERSION  ROUTINE 


350  ElEM********** *************** 
360  PD-1 

370  FOR  L-l  TO  N 
380  D-0 

390  FOR  P-1  TO  N 

400  D-D*G#(L,P)*G#(L,P) 

410  NEXT  P 
420  D-SQR(D) 

430  PD»PD*D 
440  NEXT  L 
450  DETM-1 
460  FOR  L-i  TO  N 


I 


470  Q-L 

480  J#(L+20)-Q 

490  NEXT  L 

500  FOR  L-l  TO  N 

510  CH3:M»L 

520  FOR  P-L  TO  N 

530  IF  (ABS(C)-ABS(G#(L,P)))>-0  THEN  560 

540  M-P 

550  OG#(L,P) 

560  NEXT  P 

570  IP  L-M  THEN  660 

580  R«J#(M*20) 

590  J#(fH-20)-J#(L+20) 

600  J#(L+20)»R 
610  FOR  P-1  TO  N 
620  S«G#(P,L) 

630  G#(P,L)-G#(P,M) 

640  G#(P,M)«S 

650  NEXT  P 

660  G#(L,L)»1  it 

670  D£TM-DETM*C 

680  FOR  M-l  TO  N 

690  G#(L,M)-G#(L,M)/C 

700  NEXT  M 

710  FOR  M-l  TO  N 

720  IF  L-M  THEN  790 

730  C-G#(M,L) 

740  IF  C-0  THEN  790 

750  G#(H,L)-0 

760  FOR  P-1  TO  N 

770  G#(M,P)-G#(M,P)-C*G#(L,P) 

780  NEXT  ? 

790  NEXT  M 
300  NEXT  L 
310  FOR  L-l  TO  N 
320  Q-L 

330  IF  J#(L+-20)«Q  THEN  950 
840  H-L 
350  M-H+l 

360  IF  J#(f*+-20)-Q  THEN  880 
370  IF  N>H  THEN  850 
380  J#(M+20)-J#(L+20) 

890  FOR  P-1  TO  N 
900  C»G#(L,P) 

910  G#(L,P)»G0(M,P) 

920  G#(M,P)-C 
930  NEXT  P 
940  J#(L+20)-Q 
950  NEXT  L 
960  OETM-ABS(DETH) 

970  DTNRM-DETM/PD 
930  PRINT"DTNRM-M ; DTNRH 
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REM** ****** ******************* ******************* ********** 

REM  MULTIPLIES  INVERTED  MATRIX  BY  H*K* IDENTITY  MATRIX 
REti*  ******************************************************* 

FOR  1-1  TO  N 
FOR  J-l  TO  N 
G#(I,J)-G#(I,J)*H*K 
NEXT  J 
NEXT  I 

PRINT SPRINT : PRINT-COLUMN  VECTOR  OF  EQUATION" 

^EM**************** ****** ********************************************* 
REM  CREATES  COLUMN  VECTOR  OF  EQUATION  TO  BE  EVALUATED,  MULTIPLIED  BY  h*R 

* 

PRINT"THE  FUNCTIONS  TO  BE  EVALUATED  ARE  OF  THE  PORM:  AX*2+BY~2+CX+DY+E 
PRINT-INPUT  A": INPUT  A 
PRINT- INPUT  B“: INPUT  B 
PRINT-INPUT  C“: INPUT  C 
PRINT" INPUT  0": INPUT  D 
PRINT-INPUT  £" : INPUT  £ 

FOR  i-*I  TO  N 

Q-I 

IF  ION/2  THEN  Y-1/(T+1)  ELSE  Y-2/(T+l) 

IF  K-N/2  THEN  X-Q/(SH)  ELSE  X«(Q-N/2)/(S*l ) 

F# ( I ) - ( A*X" 2+B* Y* 2+C*X*D* Y+E ) *H*K 
PRINT  F#( I) 

NEXT  I 

PRINT: PRINT: PRINT-SOLUTION  MATRIX" 

FOR  1-1  TO  N 
U#(I)-0 
FOR  J-l  TO  N 

U#(I)-U#(I)+G#(I,J)*F#(J) 

NEXT  J 

LPRINT  U#(I) 

NEXT  I 
END 


SAL ERKIN  ROUTINE  -  (ONE-DIMENSION) 


10 

20 

30 

40 

50 

60 

70 

30 

90 

100 

110 

120 

130 

140 

150 

160 

170 

130 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

430 

490 

500 

510 

520 


REM 

REM  THIS  PROGRAM  USES  A  DIRECT  GAUSSIAN  ELIMINATION  ROUTINE  tfIM 
REM  PIVOTING  TAKEN  FROM  APPLIED  NUMERICAL  METHODS  FOR  DIGITAL  COMPUTATION 
REM  Bf  JAMES,  SMITH,  &  WOLFORD,  PG  192-193.  IN  ADDITION,  THE  INITIAL 
REM  COEFFICIENT  MATRIX  HAS  BEEN  MODIFIED  SO  THAT  THE  LARGER  NUMBERS  ARE 
REM  IN  THE  UPPER  LEFT  OF  THE  MATRIX;  THIS  WAS  DONE  TO  MAXIMIZE  THE 
REM  EFFECTIVENESS  OF  THE  PIVOTING. 

REM 

PRINT  CHR$(26)+CHR$(27)+CHR$(13) 

PRINT"GAL ERKIN  (METHOD  OF  MOMENTS)  ROUTINE- 
PRINT :  PRINT :  PRINT 


PRINT-INPUT  N,  THE  NUMBER  OF  DIVISIONS"; 
INPUT  N 
LPRINT-N-  ";N 
M-N+l 


L-N-l 

DIM  A//(N,M),X/)(N),ALPHA#(N),U#(N) 
DEFD3L  Q-T,7 


REM******THIS  CALCULATES  THE  INITIAL  L  MATRIX*****4******* 
FOR  1-1  TO  N 
FOR  J-l  TO  N 


Q-M-I-.R-M-J 
A#(I,J)-(Q*R)/(Q+R+1#) 
NEXT  J 


NEXT  I 

PRINT: PRINT: PRINT-PICK  A  FUNCTION  TO  EVALUATE" 
PRINT:PRINT"1-  G-10" 

PRINT" 2-  G-X~2" 

PRINT" 3-  G-X~2  +•  1" 

PRINT"4-  G-X"2  +  X  +  1" 

PRINT" 5-  G-4X~2  +  1" 

INPUT  Z 

FOR  1-1  TO  N 

5*M-I 

ON  Z  GOSUB  960,980,1000,1020,1040 

NEXT  I 

FOR  K-l  TO  L 

H-K 

R-ABS(A#(K,K)) 

B-K+l 

FOR  I-B  TO  N 
S-ABS(A0(I,K)) 

IF  (R-S)>-0  THEN  460 

R-S:H-I 

NEXT  I 

IF  (H-K)-O  THEN  530 
FOR  J-K  TO  M 
Q«A0(H,J) 

A#(H, J )*A^(K, J ) 

A#(K,J)-Q 
NEXT  J 


530  FOR  I-B  TO  N 
540  Q-A#(I,K)/A#((C,K) 

550  FOR  J-B  TO  M 

560  A# ( I , J ) -A# ( I , J ) -Q*  A# ( K , J ) 

570  NEXT  J 
530  NEXT  I 
590  FOa  I-B  TO  M 
600  A#(l,K)-0 
610  NEXT  I 
620  NEXT  X 

630  X#(N)-A#(N,M)/A#(N,N) 

640  FOR  0-1  TO  L 

650  r-o 

660  I-N-0 

670  C-I+l 

630  FOR  J-C  TO  N 

690  T-T+A#(I,J)*X#(J) 

700  NEXT  J 

710  X#(I)-(A#(I,M)-T)/A#(I,I) 

720  MEXT  0 

730  FOR  1-1  TO  M 

740  ALPHA#(I)-X#(M-I) 

750  MEXT  I 

760  PRINT : PRINT :PRINT"INPUT  THE  POINT  YOU  WANT  EVALUATED"; 

770  INPUT  Y 
730  X-Y 

790  FOR  1-1  TO  N 
300  U#(I)-Y-Y-(I+1) 

810  NEXT  I 
320  S-0 

330  FOR  1-1  TO  M 
340  Q— U #( I )* ALPHA# ( I ) 

350  S-S+Q 
860  NEXT  I 

370  LPRIMT"U ( " ; X ; " ) -  ";S 

880  PRINT: PRINT" DO  YOU  WANT  TO  EVALUATE  ANOTHER  POINT"; 

890  INPUT  Y$ 

900  IF  Y$»"Y"  THEN  760 
910  END 

930  REM  THESE  SUBROUTINES  CALCULATE  IHE  INHOM.  TERM  TO  ADO  TO  THE  AUGMENTED 
940  REM  COEFFICIENT  MATRIX 

950  ZZlt***™****  *************  *********************  kirk*  k-trJrk*-k**k****kk*k*k 

960  A#(I,N+l)-5#*3/(S+2#) 

970  RE'TURN 

980  A#(I,N+l)-S/(4#*(S+4#)) 

990  RETURN 

1000  A#(I,M+l)-S*(3#*S+10#)/(4#*(S+4#)*(S+2#)) 

1010  RETURN 

1020  A#(I,N+l)»U3#*(S+4#)*(S+3#)*(S+2#)-12#*(3#*S"2+13#*S+26#))/(12#*(S+4#) 
*(S+3#)*(S+2#)) 

1030  RETURN 

1040  A#(I,N+l)-SM3MS+3#)/(2#*(S+2#)*(SM#)) 

1050  RETURN 
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GAL ERKIN  ROUTINE  -  (TWO-DIMENSION) 

10  REM  *********************************************************  *************** 
20  REM  TRIS  PROGRAM  USES  A  DIRECT  GAUSSIAN  ELIMINATION  ROUTINE  WITH 
30  REM  PIVOTING  TAKEN  FROM  APPLIED  NUMERICAL  METHODS  FOR  DIGITAL  COMPUTATION 
40  REM  BY  JAMES,  SMITH,  &  WOLFORD,  PG  192-193. 

50  REM  *********************************************************************** 
SO  PRINT  CHR$(26)+CHR$(27)+CHR$(13) 

70  PRINT"GALERKIN  (METHOD  OF  MOMENTS)  ROUTINE" 

30  PRINT: PRINT: PRINT 

90  PRINT" INPUT  N,  THE  NUMBER  OF  DIVISIONS"; 

100  INPUT  N 
110  LPRINT"N-  ";N 
120  M-NH 
130  L-N-l 

140  DIM  A#(N,M) ,X#(N) , ALPHA# (N) ,U0(N) 

150  DEFDBL  A-H,Q-T,Y,X 

160  REM***** ******* *******************  ************************* *********** 

170  REM******THIS  CALCULATES  THE  INITIAL  L  MATRIX***********-** 

130  FOR  1-1  TO  N 

190  FOR  J»1  TO  N 

200  Q-I:R»J:S-2*(Q+R) 

210  P-I  MOD  3:V-J  MOD  3 

220  IF  (P-1  AND  V-l)  THEN  GOSUB  1200 

230  IF  (P-1  AND  V-2)  THEN  GOSUB  1250 

240  IF  (P-2  AND  V-l)  THEN  250  ELSE  270 

250  Q-J:R-I 

260  GOSUB  1250 

270  IF  (P-1  AND  V-0)  THEN  GOSUB  1320 
230  IF  (P-0  AND  V-l)  THEN  290  ELSE  310 
290  q-J:R-I 
300  GOSUB  1320 

310  IF  (P-2  AND  V-2)  THEN  GOSUB  1390 

320  IF  (P-2  AND  V-O)  THEN  GOSUB  1440 

330  IF  (P-0  AND  V-2)  THEN  340  ELSE  360 
340  Q-J :R-I 
350  GOSU3  1440 

360  IF  (P-0  AND  V-0)  THEN  GOSUB  1510 

370  NEXT  J 

330  NEXT  I 

390  PRINT  CHR$(7) 

400  PRINT: PRINT" THE  FUNCTIONS  TO  EVALUATE  ARE  OF  THE  FORM:  AX~2+BY~2+CX+DY+£" 
410  PRINT" INPUT  A"; 

420  INPUT  A 

430  PRINT" INPUT  B"; 

440  INPUT  B 

450  PRINT" INPUT  C"; 

460  INPUT  C 

470  PRINT"INPrJT  D"; 

480  INPUT  D 

490  PRINT"INPUT  E"; 

500  INPUT  E 
510  FOR  1-1  TO  N 


530  P-Q  MOD  3 

540  IP  P-1  THEN  GOSUB  1560 
550  IP  P-2  THEN  GOSUB  1620 
560  I?  P-0  THEN  GOSUB  1680 
570  NEXT  I 
580  FOa  S-l  TO  L 
590  H-K 

600  R-ABS(A#(K,iO) 

610  B-K+l 

620  FOR  I-B  TO  N 

630  S-ABS(A#(I,K)) 

640  IP  (R-S)>-0  THEN  660 
650  R-S:H-I 
660  NEXT  I 

670  IP  (H-K)-O  THEN  730 
680  FOR  J-K  TO  M 
690  Q— A#(H,J) 

700  A#(H,J)-A#(R,J) 

710  A#(K,J)-Q 

720  NEXT  J 

730  FOR  I-B  TO  N 

740  Q-A#(I,K)/A#(K,K) 

750  FOR  J-B  TO  M 

760  A#(I,J)«A#(I,J)-Q*A#(K,J) 

770  NEXT  J 
780  NEXT  I 
790  FOR  I-B  TO  N 
800  A#(I,K)-0 
310  NEXT  I 
820  NEXT  K 

830  X#(N)-A#(N,M)/A#(N,N) 

340  FOR  0-1  TO  L 

350  T-0 

860  I-N-0 

370  C-I+l 

830  FOR  J-G  TO  N 

890  T*T+A#(1,J)*X#(J) 

900  NEXT  J 

910  X#(I)-(A#(I,M)-T)/A#(I,I) 

920  NEXT  0 

930  FOR  1-1  TO  N 

940  ALPHA#(I)-X#(I) 

950  NEXT  I 

960  PRINT  CHR$ ( 7 ) 

970  PRINT: PRINT: PRINT" INPUT  THE  X,Y  POINTS  YOU  WANT  EVALUATED 

980  INPUT  X,Y 

990  FOR  1-1  TO  N 

1000  Q-I 

1010  P-I  MOD  3 

1020  IF  P-1  THEN  U#(I)-(X*Y)"((2/if*q+l#)/3#)*(l#-X)*(l#-Y) 

1030  IF  P-2  THEN  U#(I)-X"((2/^+5#)/3#)*Y*(l#-X)*(l#-Y) 

1040  IF  P-0  THEN  U#(I)-Y"((2#*Q+3#)/3«»)*X*(l//-X)*(l#-Y) 

1050  NEXT  I 
1060  S-0 


1070  FOR  1-1  TO  N 
1080  Q-U#(I)*ALPHA#(I) 

1090  S-S+Q 
1130  NEXT  I 

1110  LPRINT  "U(“;X;Y;")-";S 

1120  PRINT: PRINT"DO  YOU  WANT  TO  EVALUATE  ANOTHER  POINT"; 

1130  INPUT  Y$ 

1140  IF  Y$-"Y"  THEN  970 
1150  END 

1160  RSM**^***^***^*^**^***-****** k*k ****** kkkkkk***kk*kk**kk*kkkkkkirk*** 

1170  REM  THE  FOLLOWING  SUBROUTINES  ARE  USED  TO  CREATE  THE  AUGMENTED  LMN 
1130  REM  COEFFICIENT  MATRIX 

1200  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-1,4,7,...  AND  N-l ,4,7, . . . 

1210  A»l*/(S+5#)-2#/(S+3#)+l#/(S+ll#> 

1220  B-(R-l#)/(S-l#)-(2#*R+l#)/(S+2#)+(R+2#)/(S+5#) 

1230  A#(I,J)-4#*(2#*R+1#)*A*8 
1240  RETURN 

1250  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-1,4,7,...  AND  M-2,5,3,... 

1260  A»l)//(2#*Q+7//)-2#/(2#*0+-10//>)-H#/(2//*Q+13#) 

1270  8-(R+1#)/(S+3#)-(2MR45#)/(S+60)+(R+40)/(S+9iP) 

1230  C»l///(S+9#)-2#/(S+12#)+l#/(S+150) 

1290  D»l#/(2#*Q+7#)-l#/(2#*Q+4#) 

1300  A#(I,J)-(2#*A*B*(2#*R+5#)+18#*C*D) 

1310  RETURN 

1320  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-1,4,7,...  AND  N-3,6,9,... 

1330  A-l#/(2#*Q+7#)-2#/(2#*Qfl0#)+l#/(2#*Q+13#) 

1340  3-R/(SM#)-(2#*a+3#)/<3+4#)+<R+3#)/(3+7#> 

1350  C-l#/(S+7#)-2#/(S+10#)+l#/(S+13#) 

1360  D»l#/(2#*Q+7#)-l#/(2#*Q+4#) 

1370  A#(I,J)-(2#*A*B*(2#*R+3#)+l3#*C*D) 

1330  RETURN 

1390  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-2,5,3,...  AND  N-2,5,8,... 

1400  A-(R+l#)/(S+7#M2#*a+5#)/(S+10#)+(R+4#)/(S+13#) 

1410  B»l#/(S-H3#)-20/(SH6#)+l4y(S+19 #) 

1420  A#(I, J)-(A*(2#*R+5#)/45#-B) 

1430  RETURN 

1440  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-2,5,8,...  AND  N-3,6,9,... 

1450  A-l^/(2#*CH-ll#)-2#/(2#*Q+14#)+l#/(2#*(^17#) 

1 460  B-R/ ( 2#*R+3#)- ( 2#*R+3# ) / ( 2#*R+6# )+( R+3# ) / ( 2#*R+9# ) 

1470  C-U/(2#*R+9#)-2#/(2#*R+l2#)+10/(2#*R+15#) 

1480  D-1#/(2#*Q+11#)-1#/ (2#*Q+30) 

1490  A#(I,J)«(2#*A*B*<2#*R+3#)+18£*C*D) 

1500  RETURN 

1510  REM  CALCULATES  LMN  MATRIX  VALUES  FOR  M-3,6,9,...  AND  N-3,o,9,... 

1520  A-R/(S+3*0-(2//*R+3#)/(S+6#)+(R+6#)/(S+9$) 

1530  B-ltf/(S+9#)-2/>/(S+12#)*l#/(S+150) 

1540  A#(I,J)-(A*(20*R+3#)/45#-8) 

1550  RETURN 

1560  REM  CALCULATES  GM  FOR  M-1,4,7,... 

1570  F-(A48)*(l#/(R+10#)-l#/(a+l3tf)) 

1530  G-<C+D)*(l#/(R+7#)-l#/(a+10#)) 

1590  H-E*(l#/((R+4/n*(a+4#))-2#/<(R+40)*(R+70))-H0/<(a+7f?)*<R+7#)>) 

1600  A#(I,N+l)-9M(F-H3)*(l#/(R+4#)-10/(R*7#))+H) 


1610  RETURN 

1620  REM  CALCULATES  GM  FOR  M»2,5,8,... 

1630  F»A/2#*(l#/(R+14#)-10/(a+17<O) 

1640  G-C/2i7*(l#/(R+110)-10/(R+14//)) 

1650  H*(l#/ (R+3#)-l#/ (R-H1#) )*(3#*B/20ijf+D/ 4//+E/2#) 
1660  A#  ( I ,  N-H  )»F+GH1 
1670  RETURN 

1680  REM  CALCULATES  GM  FOR  M-3,6,9,... 

1690  ?»(l#/(a+6#)-l#/(a+9i?))*(3#*A/20#-HC/4#+£/2#) 
1700  G-B/2#*(1#/(R+12#)-1#/(R-H5#)) 

1710  ri-D/2#*(l#/(R+9#)-l#/(R+12#)) 

1720  A#(I,N+1)-F+GH1 
1730  RETURN 


COLLOCATION  ROUTINE  -  (ONE-DIMENSION) 


10  ASM  ****************************** *** ************ ************** ******** 
20  REM  THIS  PROGRAM  USES  AN  INVERSION  ROUTINE  TAREN  PROM  NUMERICAL  METHODS 
30  REM  ST  R.W.  HORN BECK,  PG  294-295.  THE  ROUTINE  EMPLOYS  GAUSS-JORDAN 
40  REM  ELIMINATION  WITH  COLUMN  SHIFTING  TO  MAXIMIZE  PIVOT  ELEMENTS. 

50  REM  *************4****** *************** ****** ************************** 
60  PRINT  CHR$(26)+CHR$(27)+CHR$(13) 

70  LPAINT"CO-LOCATION  (METHOD  OF  MOMENTS)  ROUTINE** 

80  PRINT: PRINT: LPRINT 

90  PRINT" INPUT  N,  THE  NUMBER  OF  DIVISIONS"; 

100  INPUT  N 
110  LPRINT" N»  ";N 

120  DIM  C#(N,N),ALPHA#(N),G#(N),J#(N+25),U#(N) 

130  DSFDBL  Q-S,Y,C,D 

150  REM  C#(N,N)  CONTAINS  THE  INITIAL  L  MATRIX,  THEN  THE  INVERTED  L  MATRIX; 
160  REM  J#(N+25)  IS  USED  ONLY  IN  THE  INVERSION  ROUTINE 

1 70  REM ******************************************************************* 
180  REM****** THIS  CALCULATES  THE  INITIAL  L  MATRIX**********-*** 

190  FOR  1-1  TO  N 
200  FOR  J-l  TO  N 
210  Q-I:R»J:S-N 

220  C#(I,J)-R*(R+l#)*(q/(S+l#))"(R-l#) 

230  NEXT  J 
240  NEXT  I 
250  PD-1 

260  FOR  L-l  TO  N 
270  D-0 

280  FOR  X-l  TO  N 

290  D-D*C#(L,K)*C#(L,X) 

300  NEXT  K 
310  D-SQR(D) 

320  PD-PD*D 
330  NEXT  L 
340  DETM-1 
350  FOR  L-l  TO  N 
360  Q-L 

370  J#(L+20)-Q 
330  NEXT  L 
390  FOR  L-l  TO  N 
400  C-0 
410  M-L 

420  FOR  R-L  TO  N 

430  IF  (ABS(C)-ABS(C0(L,K)))  >-  0  THEN  460 

440  M-K 

450  C-C#(L,K) 

460  NEXT  R 

470  IF  L-M  THEN  560 

430  a-J#(M*20) 

490  J#(M+20)-J#(L+20) 

500  J#(L+20)»R 
510  FOR  X»i  TO  N 
520  S-C#(K,L) 


530  C#(R,L)-C#(X,M) 

540  C#(K,M)-S 
550  NEXT  K 
560  C#(L,L)-1# 

570  DETM«DETM*C 
530  FOR  M-l  TO  cl 
590  C#(L,M)»C#(L,M)/C 
600  NEXT  M 
610  FOR  M-l  TO  N 
620  IF  L-N  THEN  690 
630  C-C#(M,L) 

640  IF  C-0#  THEN  690 

650  C#(N,L)-0 

660  FOR  R»1  TO  N 

670  C#(M,X)«C#(M,K)-C*C#(L,X) 

680  NEXT  K 
690  NEXT  M 
700  NEXT  L 
710  FOR  L»1  TO  N 
720  Q-L 

730  IF  J#(L+-20)-Q  THEN  850 
740  M-L 
750  M-M+l 

760  IF  J#(M+20)-Q  THEN  730 
770  IF  N>M  THEN  750 
780  J#(M+,20)-J#(L+20) 

790  FOR  E-l  TO  N 
300  C«C#(L,X) 

310  C#(L,K)-C#(M,K) 

820  C#(M,K)-C 
330  NEXT  K 
340  J#(L+20)-Q 
850  NEXT  L 
360  DETM-A8S ( DETM ) 

370  DTNRM-DETM/PD 

330  REM* ********  CALCULATES  gm,  THE  INNER  (PRODUCT  OF  THE  WEIGHT  FN  &  g ******* 
390  FOR  1-1  TO  N 
900  Q-I:R-N 

910  REM*****THE  NEXT  LINE  CAN  BE  CHANGED  FOR  OTHER  INHOMOGENEITIES****** 

920  3#(I)-l#+4#*(Q/(a+l#))~2 
930  PRINT"G(";I;")-  ";G#(I) 

940  NEXT  I 

950  REM******CALCULATES  ALPHA******** 

960  FOR  I»l  TO  ft 
970  ALPHA# ( I) -0 
930  FOR  J-l  TO  N 

990  ALPHA#(I)-ALPHA#(I)+C#(I,J)*G#(J) 

1000  NEXT  J 

1010  PRINT" ALPHA ( " ; I; " )»  ";ALPHA#(I) 

1020  NEXT  I 

1030  PRINT" INPUT  THE  THE  POINT  YOU  WANT  EVALUATED"; 

1040  INPUT  Y 
1050  X-Y 

1060  FOR  1-1  TO  N 
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1370  U#( I)*Y-Y"( 1+1 ) 

1030  PRINT  U0(I) 

1090  NEXT  I 
1100  S-0 

1110  FOR  1-1  ro  n 
1120  Q*U0(I)*ALPHA#(I) 

1130  S-S+U#(I)*ALPHA#(I) 

1140  PRINT"U(“;I;")»  ";Q 
1150  NEXT  I 

1160  PRINT"U(";X;")-  ";S 
1170  LPRINT"U ( " ; X; " )»  ";S 

1180  PRINTS  PRINT" DO  YOU  WANT  TO  EVALUATE  ANOTHER  POINT  (Y,N)" 
1190  INPUT  Y$ 

1200  IP  Y$-"Y"  THEN  1030 
1210  END 


COLLOCATION  ROUTINE  -  (TWO-DIMENSION) 


10  REM  ********************************************************************* 

20  REM  THIS  PROGRAM  USES  A  DIRECT  GAUSSIAN  ELIMINATION  ROUTINE  WITH 
30  REM  PIVOTING  TAKEN  FROM  APPLIED  NUMERICAL  METHODS  FOR  DIGITAL 
40  REM  COMPUTATION  BY  JAMES,  SMITH  &  WOLFORD,  PG  192-193. 

60  PRINT  CHR$ ( 26 ) +CHR$ (27 )+CHR$ (13) 

70  PRINT: PRINT: PRINT 
30  DEFDBL  T,Z 

90  PRINT~IN?UT  THE  NUMBER  OF  X  DIVISIONS,  AND  THE  NUMBER  OF  Y  DIVISIONS" 

100  INPUT  T,Z 
110  N»Z*T 
120  M-N+1:L-N-1 
130  LPRINT"N»";N 

140  DIM  C#(N,M),AL?HA#(N),X#(N),U#(N) 

150  DEFDBL  Q-S,X,Y ,A-£, V,T 

160  REM******************************************************************* 
170  REM****** THIS  CALCULATES  WE  INITIAL  L  MATRIX***** ********  *********** 

210  FOR  1-1  TO  M 
220  V-I 

230  IF  K-N/2  THEN  Y-1/(Z+1)  ELSE  Y-2 /(Z+l) 

250  IF  K-N/2  THEN  X-V/(T+1)  ELSE  X«( V-N/2)/(T+l ) 

270  FOR  J-l  TO  N 

280  Q-J 

290  P-J  MOD  3 

300  IF  ?-l  THEN  GOSUB  1140 
310  IF  P-2  THEN  GOSUB  1210 
320  IF  P-0  THEN  GOSUB  1260 
340  NEXT  J 
360  NEXT  I 

370  REM********* CALCULATES  gra,  THE  INNER  PRODUCT  OF  THE  WEIGHT  FN  &  g ******* 
375  PRINT  CHR$(7) 

330  PRINT :?RINT"THE  FUNCTIONS  TO  BE  EVALUATED  ARE  OF  WE  FORM:  AX~2  +  BY~2 
+  CX  +  DY  -HE" 

390  PRINT" INPUT  A"; 

400  INPUT  A 

410  PRINT" INPUT  B"; 

420  INPUT  3 

430  PRINT"INPUT  C"; 

440  INPUT  C 

450  PRINT" INPUT  0"; 

460  INPUT  D 

470  PRINT" INPUT  E"; 

480  INPUT  E 
490  FOR  1-1  TO  N 
500  Q-I 

510  IF  K-N/2  THEN  Y-1/(Z+1)  ELSE  Y-2/(Z+l) 

520  IF  K-N/2  THEN  X-Q/(T+1)  ELSE  X-(Q-N/2)/( T+l ) 

530  C#(I,M)-A*X*X  *  B*T*Y  +  C*X  +  D*Y  +  E 

540  NEXT  I 

550  FOR  K-l  TO  L 

560  H-K 
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570  R-ABS(C#(X,lD) 

580  B-X+l 

590  FOR  I-B  TO  N 

600  S-ABS(C#<I,K)) 

610  IF  (R-S)>-0  THEN  630 
620  R-S:H-I 
630  NEXT  I 

640  IF  (H-K)-O  THEN  700 
650  FOR  J-K  TO  M 
660  Q«C#(H,J) 

670  C#(a,J)-G#(K,J) 

680  C#(K,J)-Q 

690  NEXT  J 

700  FOR  I-B  TO  N 

710  Q«C#(I,K)/C#(R,X) 

720  FOR  J-B  TO  M 

730  C#(I,J)-C#(I,J)-Q*C#(K,J) 

740  NEXT  J 
750  NEXT  I 
760  FOR  I-B  TO  N 
770  C#(I,K)-0 
780  NEXT  I 
790  NEXT  K 

800  X#(N)-C0(N,M)/C#(N,N) 

910  FOR  0-1  TO  L 

320  S-0 

830  I-N-0 

340  C-I+l 

350  FOR  J-G  TO  N 

360  S»S+C#(I,J)*X#(J) 

370  NEXT  J 

380  X#(l)-<C#(I,M)-S)/C#(I,t) 

390  NEXT  0 

900  FOR  I-i  TO  N 

910  ALPrfA#(I)-X#(I) 

920  NEXT  I 

925  PRINT  CHR$ ( 7 ) 

930  PRINT" INPUT  THE  X,Y  POINT  YOU  WANT  EVALUATED"; 

940  INPUT  X,Y 
950  FOR  1-1  TO  N 
960  q-i 
970  P-I  NOD  3 

980  IF  P-1  THEN  U#(I)-(X*Y)*((2#*Q+1#)/3#)*(1#-X)*(U-Y) 
990  IF  P-2  THEN  U#(l)-X~((2#*Q+5#)/3#)*Y*(l#“X>*(l#-Y> 
1000  IF  P-0  THEN  U#(I)-Y*((2#*Q+3#)/3#)*X*(l#-X)*(l#-Y) 
1010  NEXT  I 
1020  S-0 

1030  FOR  1-1  TO  N 
1040  Q-U#(I)*ALPHA#(I) 

1050  S-S+UiK I)*ALPHA# (I) 

1060  ?RINT"U(";I;")»  ";Q 
1070  NEXT  I 

1080  PRINT"U(";X;Y;")«  ";S 
1090  LPRINT"U(";X; Y;")«  ";S 


lioo  print:print"do  you  want  to  evaluate  another  point  (Y,N) 
1110  INP'JT  Y$ 

1120  IF  YI-^'Y"  THEN  930 
1130  END 

1140  REM  CALCULATES  THE  LMN  MATRIX  VALUES  FOR  M-1,4,7,... 
1150  A“X'*((2#*Q-5#)/3#)*(Y‘*((20*Q+l#)/3#)-Y“((2#*Q+4#)/3#) ) 
1160  B-Y*((2#*Q-5#)/3#)*(X-((2#*Q+l#)/3#)-X^((2#*Q+4#)/3#)) 
1170  C-X*((2#*Q-2#)/3#)*(Y~((2#*Q+4#)/3#)-Y“((2#*Q+l#)/3#)) 
1130  OY~((2#*Q-2#)/3#)*(X-((2#*Q+4#)/3#)-X~((20*Q+l#)/3#)) 
1190  C#(I,J)-((2#*Q-2#)*(A+B)+(2#*Q+4#)*(C+-D))*(2#*Q+l#)/9# 
1200  RETURN 

1210  REM  CALCULATES  THE  LMN  MATRIX  VALUES  FOR  N-2,5,3,... 
1220  A»(2#*Q+2#)*X-((2#*Q-l#)/3#)-(2#*Q+8#)*X“((2#*Q+2#)/3#) 
1230  B*2#*(X~((2#*Q+8#)/3#)-X‘‘((2#*Q+5#)/3#)  ) 

1240  C#(I,J)»A*<2#*Q+5#)*(Y-Y*7)/9#+3 
1250  RETURN 

1260  REM  CALCULATES  THE  LMN  MATRIX  VALUES  FOR  N-3,6,9,... 
1270  A-(2#*Q)*Y*((2#*Q-3#)/3iif)-(2#*Q+6#)*Y*(2iif*Q/3#) 

1230  B-2#*(Y-((2#*Q+6#)/3^>-r((2#*Q+3#)/3#)) 

1290  C//(I,J)«A*(2#*Q+3#)*(X-X*X)/9#+B 
1300  RETURN 
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The  purpose  of  this  study  was  to  determine  the  feasi¬ 
bility  of  using  the  method  of  weighted  residuals  to  obtain 
approximations  to  the  discrete  Green's  function,  or  analogs 
to  it.  The  weighted  residual  methods  of  Galerkin  and  collo¬ 
cation,  ad  well  as  the  finite  difference  method  were  programed 
on  a  Kaypro  II  micro-computer  in  Microsoft  Basic.  These 
programs  were  used  to  generate  approximations  to  the  one- 
and  two-dimensional  Poisson's  equation.  The  two-dimensional 
case  was  restricted  to  the  geometry  of  a  unit  square. 

Various  inhomogeneity  terms  were  used  to  obtain  approximate 
solutions  to  the  discrete  Green's  functions  or  their  analogs. 
The  results  were  compared  with  the  analytical  values  at 
various  interior  nodal  points  on  the  mesh;  The  average  percent 
error  for  the  approximations  were  reported  for  each  case  as 
the  number  of  interior  nodal  points  was  increased.  The  areas 
of  consideration  were  the  rate  of  convergence  to  the  analytical 
solution,  the  amount  of  time  it  took  to  run  each  program,  and 
the  accuracy  of  the  approximate  solutions.  The  results  of 
this  study  indicate  that  the  Green's  functions  or  analogs 
obtained  are  valid  approximations  to  the  discrete  Green’s 
functions.  The  method  of  weighted  residuals  proved  to  be 
very  sensitive  to  the  choice  of  basis  functions,  resulting 
in  ill-conditioned  matricies  in  some  instances. 
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